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ABSTRACT 


This  report  describes  31  investigations  in  the  field  of  seismic 
discrimination.  These  are  grouped  as  follows:  body- wave 

studies  (6  contributions),  surface- wave  studies  (5  contribu- 
tions), studies  of  long-  and  short-period  magnitude,  explosion 
yield,  and  seismic  energy  (9  contributions),  miscellaneous 
studies  (9  contributions),  and  recent  developments  in  our  data 
and  computer  systems  (2  contributions). 
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SUMMARY 


This  is  the  twenty-sixth  Semiannual  Technical  Summary  report  describing  the  activities  of 
Lincoln  Laboratory,  M.I.T.,  in  the  field  of  seismic  discrimination.  Due  to  the  recent  change 
in  fiscal  year,  this  report  covers  the  period  1 July  1976  to  31  March  1977.  The  objective*  of  the 
Lincoln  Laboratory  program  is  to  carry  out  fundamental  research  into  the  seismological  prob- 
lems associated  with  the  detection,  location,  and  identification  of  earthquakes  and  nuclear  ex- 
plosions, and  with  the  estimation  of  explosion  yields.  In  order  to  investigate  these  probltems, 
we  are  continuously  improving  our  ability  to  manipulate  and  process  seismic  data  from  a global 
network  of  high-quality  digital  stations  and  arrays. 

During  the  period  covered  by  this  report,  data  have  become  available  from  the  firs  t rSeismic 
Research  Observatories  (SROs).  At  the  time  of  writing,  data  are  being  received  on  a routine 
basis  from  six  of  these  stations:  ANMO  (Albuquerque),  MAIO  (Mashad),  GUMO  (Guam),  NWAO 
(Australia),  SNZO  (New  Zealand),  and  TATO  (Taiwan).  Fairly  extensive  data  are  now  available 
for  the  first  three  of  these,  and  several  contributions  in  this  report  describe  our  first  attempts 
to  use  these  new  data  for  seismological  research.  The  quality  of  these  data  is  very  high,  and 
we  expect  them  to  make  a major  impact  on  seismology  as  more  stations  become  operational. 

In  the  area  of  general  body-wave  studies,  several  investigations  are  described  in  this  re- 
port. We  note  the  presence  of  small  emergent  phases  ahead  of  larger,  more-impulsive  phases 
in  many  cases  of  shallow  earthquakes,  using  the  short-period  (SP)  SRO  data.  This  cleatly 
complicates  the  selection  of  onset  times  from  conventional  WWSSN  film-chip  data.  Other  studies 
have  used  P-wave  travel  times  to  investigate  the  variation  in  crust  and  upper-mantle  structure 
across  the  United  States,  and  to  determine  velocities  in  the  source  region  by  a master  event 
technique.  It  appears  possible  to  directly  measure  P-  and  S-wave  velocity  in  source  ar£as  such 
as  downgoing  lithosphere  slabs,  and  make  conclusions  about  anomalous  zones  within  thepe  slabs. 
Earlier  work  on  the  representation  of  the  seismic  source  by  the  moment  tensor  has  been  ex- 
tended to  body- wave  sources. 

Surface- wave  studies  described  include  a new  rapid  algorithm  for  the  determination  of  group 
and  phase  velocity,  a series  of  observations  of  Rayleigh-wave  dispersion  at  the  Mashad  SRO, 
and  a determination  of  structure  in  the  Novaya  Zemlya  area  using  a pair  of  events.  In  another 
study,  observations  of  NTS  explosions  at  Albuquerque  SRO  are  analyzed  for  group  velocity  dis- 
persion using  a data  adaptive  prediction  error  algorithm.  A study  of  surface-wave  overtones  is 
also  included,  and  a technique  is  described  for  the  problem  of  mode  separation. 

A series  of  investigations  continuing  our  research  into  the  determination  of  magnitude  and 
yield  are  described.  A study  of  the  variation  with  azimuth  of  the  long-period  (LP)  P-wfcve 
amplitudes  from  two  explosions  shows  features  that  are  consistent  with  an  earlier  study  of  sta- 
tion magnitude  bias  using  short-period  data.  An  initial  application  of  SRO  data  to  the  pr  oblem 
of  seismic  scaling,  and  a new  technique  for  measuring  the  amplitude  and  phase  response  of  the 
WWSSN  SP  instrument  are  described.  We  have  continued  our  investigation  of  methods  for  the 
determination  of  station-detection  characteristics.  We  are  particularly  concerned  about  the 
inclusion  of  saturation  parameters,  and  the  proper  evaluation  of  the  parameters  describing 
seismicity.  A series  of  stations  with  known  parameters  are  then  used  as  a model  of  a global 
network,  and  a synthetic  earthquake  catalog  for  the  network  has  been  generated  using  a simula- 
tion program.  A theoretical  development  describes  how  the  seismic  energy  radiated  by  a 
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finite  propagating  line  source  can  be  evaluated  using  the  moment  tensor  formulism,  and  a sim- 
ilar approach  is  used  to  calculate  the  explosive  moments  and  radiated  Rayleigh- wave  energies 
for  a set  of  U.S.  explosions.  In  another  study,  the  LP  P- waves  observed  from  Longs  hot,  Milrow, 
and  Cannikin  are  compared  with  observed  SP  data.  It  appears  that  LP  data  may  provide  a much 
more  stable  measure  of  yield. 

A wide  variety  of  miscellaneous  studies  is  included.  These  range  from  an  extension  of 
maximum  entropy  spectral  analysis  to  the  study  of  multichannel  cross  spectra,  to  studies  of  Q 
structure  in  the  earth.  Cepstral  analysis  is  applied  to  an  event  near  Semip  ala  tins  k,  and  it  is 
shown  to  be  a multiple  event.  An  attempt  to  use  the  same  analysis  to  detect  the  splitting  of  the 
pP  phase  due  to  double  reflection  in  oceanic  areas  is  described.  The  application  of  polarization 
filtering  to  LP  SRO  data  is  discussed,  and  it  appears  to  be  a promising  technique  for  separating 
interfering  phases.  A theoretical  development  for  the  excitation  of  normal  modes  by  a propa- 
gating fault  is  included.  The  determination  of  the  dissipation  parameter  Q is  of  considerable 
current  interest.  Direct  measurement  of  Q for  a number  of  normal  modes  of  the  earth  is  out- 
lined, and  the  effects  of  a frequency-dependent  Q on  velocity  dispersion  are  discussed.  An 
attempt  at  carrying  out  a simultaneous  inversion  of  normal-mode  data  to  obtain  both  elastic  and 
inelastic  parameters  is  described,  using  finite  strain  theory.  A possible  correlation  between 
variations  in  the  level  of  seismic  activity  with  time  and  changes  in  the  rate  of  rotation  of  the 
earth  is  discussed. 

We  continue  to  develop  our  data-handling  facilities.  Our  PDP-11  system  is  now  a general- 
purpose  time-sharing  facility,  and  we  are  building  up  our  applications  software  library  for 
seismic  research.  The  system  allows  us  to  connect  with  the  ARPANET,  to  larger  computers, 
and  to  the  datacomputer.  Interaction  with  the  datacomputer  using  newly  developed  software  is 
increasing. 


M.  A.  Chinnery 
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I.  BODY-WAVE  STUDIES 

A.  SHORT-PERIOD  P-WAVEFORMS  FROM  LARGE  EVENTS 

RECORDED  AT  SRO  SITES 

As  part  of  an  evaluation  of  the  SRO  system,  a study  has  been  made  of  short-period 
P-waveforms  recorded  at  SRO  sites  from  large  (PDE  m^  > 5.6)  events  during  January  a 
February  1976.  The  stations  operating  at  this  time  were  AN  MO  (Albuquerque),  MAIO  (M^shad), 
and  GUMO  (Guam).  The  data  quality  at  ANMO  and  MAIO  is  generally  excellent,  but  at  QUMO 
extremely  high  noise  levels  (predominant  period  ~1  sec,  amplitude  frequently  exceeding  100  m}i) 
severely  degrade  recording  capability.  Noise  levels  at  ANMO  and  MAIO  appear  to  be  generally 
less  than  2 mp  in  the  SP  band. 

This  difference  in  noise  levels  is  reflected  in  the  detection  capability  of  the  stations-  The 
SP  data  are  saved  only  when  triggered  by  the  detection  algorithm.  There  were  55  events  of 
PDE  m^  >5.6  during  the  first  2 months  of  1976:  restricting  observations  to  only  events  of 

20  < A < 100°,  the  number  of  these  which  should  have  been  detected  at  ANMO,  MAIO,  and 
GUMO  were  2 0,  19,  and  33,  respectively.  Actual  detections  were  19,  17,  and  9,  discarding 

GUMO). 
at 


periods  during  which  the  stations  were  not  operating  (3  days  at  ANMO,  6 at  MAIO,  0 at 
The  above  give  detection  rates  for  events  of  m^  > 5.6  of  95  percent  at  ANMO,  89  percen 
MAIO,  and  a distressingly  low  27  percent  at  GUMO. 

The  nature  of  the  waveforms  recorded  is  also  of  interest.  Figure  1-1  shows  P-wav^s  [(1)  — 
(20)]  recorded  at  the  SRO  sites  from  sources  of  varying  source  depth,  distance,  and  siz£. 
Waveforms  from  deep  (h  > 70  km)  earthquakes  [(1)  — (7)]  all  have  extremely  sharp  and  impulsive 
onsets.  Those  from  shallow  events  often  show  a small  emergent  phase  up  to  4 sec  prior*  to  a 
more  impulsive  larger  arrival  [(8),  (12),  (14),  (15),  and  (19)].  The  larger  arrival  is  probably 
the  stopping  phase,  which  occurs  virtually  simultaneously  on  the  edges  of  the  fault  plane;  the 
smaller  precursive  arrival  is  due  to  the  onset  of  faulting  at  a point  of  nucleation. 

Previous  work  at  Lincoln  Laboratory  attempted  to  predict  such  complicated  waveforms 

2 

on  WWSSN  records  using  a Savage-type  source;  the  advent  of  the  SROs  gives  us  a high-rquality 
digital  dataset  which  we  intend  to  test  against  the  source  model.  An  interesting  feature  of  these 
complicated  waveforms  is  that  the  gradual  onset  may  frequently  be  masked  by  noise,  artd  that 
time  picks  made  from  such  data  may  consequently  often  be  in  error  by  several  seconds 

R.G.  North 


B.  P-WAVE  TRAVEL-TIME  ANOMALIES  AND  GRAVITY  ANOMALIES  IN  A STUDY 

OF  THE  UPPER  MANTLE  BENEATH  THE  UNITED  STATES 

3-5 

According  to  recent  studies,  gravity  anomalies  in  the  wavelength  range  of  500  to  2000  km 
may  well  be  due  to  lateral  inhomogeneities  in  density,  to  depths  of  a few  hundred  kilometers. 

On  the  other  hand,  the  analysis  of  station  anomalies  for  body-wave  travel  times  on  a sifnilar 
scale^'7  shows  that  the  corresponding  velocity  anomalies  must  originate  deeper  than  in  the 
crust. 

This  study  is  part  of  a more  extensive  project  attempting  to  correlate  gravity  data  and  body- 


wave  travel-time  data  in  relation  to  the  structure  of  the  upper  mantle  beneath  the  continental 
United  States. 


In  order  to  invert  travel  times  for  P-velocity  structure,  the  upper  mantle  of  the  U.  S.f  to 
a depth  of  800  km,  was  divided  into  5°  X 5°  blocks  of  adjustable  thickness,  and  a selection  of 
P-wave  travel-time  data  was  made  so  as  to  optimize  the  sampling  of  these  blocks  by  the  corre- 
sponding rays.  From  the  ISC  catalogs  for  the  years  1964  to  1973,  102  seismic  stations  were 
selected  in  the  U.S.  and  southern  Canada,  on  the  basis  that  they  reported  more  than  200  events 
in  that  period  of  time.  Then  2653  events  were  chosen,  using  the  following  criteria: 

Magnitude  m^  ^5.8,  event  observed  by  at  least  50  stations  in  the  world 

Event  observed  by  at  least  15  North  American  stations  (in  order  to  re- 
move the  source  effect,  common  to  all  stations),  in  the  distance  range 
of  30*  to  100° 

Accurate  location  of  epicenter 

standard  error  in  latitude  < 0.05° 
standard  error  in  longitude  ^0.06° 

No  further  constraints  on  focal  depth  or  origin  time  determination,  as 
these  effects  are  common  to  all  stations. 

The  travel-time  anomalies  were  then  computed,  as  referred  to  two  different  laterally  homo- 

Q o 

geneous  models:  Jeffreys  -Bullen,  and  "PEM"  of  Dziewonski,  Hales,  and  Lap  wood.  Anomalies 
of  magnitude  greater  than  ±2  sec  were  considered  as  bad  readings  and  rejected.  Then,  the 
station  anomalies  were  computed  by  taking  out  the  average  over  all  stations  for  each  particular 
event. 

For  each  station,  an  analysis  was  made  of  the  dependence  of  the  anomaly  on  azimuth  from 
station  and  on  distance  from  station  to  event.  For  that  purpose,  the  anomalies  were  averaged 

7 

over  azimuth  ranges  of  45°  and  distance  ranges  of  10°.  As  expected  from  previous  studies, 
three  major  regions  can  be  distinguished  in  the  U.S.  [see  Figs.  1-2 (a)  and  (b)]:  strongly  positive 
anomalies  in  the  west,  strongly  negative  in  the  center,  and  a more  complicated  pattern  in  the 
east. 

By  drawing  profiles  across  aligned  stations  and  comparing  the  azimuthal  dependence  be- 
tween stations  on  a profile,  one  can  attempt  to  put  some  bounds  on  the  depths  and  locations  of 
the  velocity  anomalies. 

NW-SE  and  W-E  profiles  across  western  and  central  U.S.  suggest  depths  as  large  as  600  to 
700  km  for  the  large  fast  anomaly  in  the  central  U.S.,  shallower  and  weaker  towards  the  east, 
and  extending  south  into  the  Gulf  of  Mexico. 

1 0 

The  next  step  in  this  study  is  to  perform  the  mathematical  inversion  of  these  travel-time 
anomalies  for  velocity  structure  and  relate  them  to  the  gravity  data  and  density  structure. 

B.A.  Romanowicz 


C.  TRAVEL-TIME  RESIDUALS  FOR  P PHASES: 

A SEARCH  FOR  SOURCE  ANOMALIES 

P and  S velocities  in  mantle  seismic  zones  are  anomalously  high  by  6 to  7 percent  when 

11  1213 

averaged  over  nearly  the  entire  length  of  the  Japanese  zone  or  the  Tonga  Zone.  ’ There 
is  a growing  need  to  know  how  this  anomaly  is  distributed  along  the  slab-like  bodies  containing 
the  seismic  earthquakes,  and  whether  or  not  the  anomaly  exists  below  the  depth  of  the  deepest 
earthquakes.  This  depth  varies  between  650  and  700  km  for  the  deepest  zones. 
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In  situ  P-wave  velocities  reported  in  Fitch  ' and  in  this  SATS  show  that  positive  anom- 
alies as  large  as  10  percent  exist  in  the  region  of  the  deepest  mantle  earthquakes.  In  contrast, 

no  significant  anomaly  was  resolved  in  the  depth  range  from  450  to  550  km  in  the  Tonga  Zone. 

15 

Fitch  suggested  that  the  deep  positive  anomaly  resulted  from  an  elevated  "650  km  discon- 
tinuity,n a major  transition  zone  that  can  be  taken  as  the  boundary  between  the  upper  and  lower 
mantles.  If  this  interpretation  were  correct,  then  the  deepest  earthquakes  would  be  below  the 
top  of  this  transition  zone  and  thus  in  the  lower  mantle.  The  question  to  be  resolved  is  Whether 
or  not  the  anomalous  material  containing  the  mantle  earthquakes  descends  to  great  depth  into 
the  lower  mantle,  albeit  without  seismic  activity. 

The  seismological  approach  to  this  problem  has  been  to  analyze  travel-time  residuals. 
Jordan  and  Lynn^  and  Engdahl*7  reported  residuals  of  ScS  and  P times,  respectively,  that 
could  be  explained  by  lithospheric  material  that  had  descended  into  the  lower  mantle.  More 
recently,  Jordan18  has  proposed  the  existence  of  a slab-like  body  extending  to  possibly  1000  km 
beneath  the  sea  of  Okhotsk.  He  analyzed  residuals  of  S,  ScS,  and  ScS-S  times  from  one  deep 
earthquake. 

We  have  begun  a study  similar  to  Jordan's;  however,  our  data  are  the  times  of  P phases 
reported  in  the  Bulletin  of  the  International  Seismological  Center  (ISC).  This  data  base  Includes 
times  for  the  years  1964  through  1973.  Because  most  deep  seismic  zones  are  nearly  linear  for 
distances  of  more  than  1000  km  along  the  strike  of  the  various  island  arcs,  residuals  from  each 
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zone  can  be  superposed  on  a single  focal  sphere  as  was  done  for  the  Tonga  earthquakes 

Figs.  1-3 (a)  through  (c).  About  40,000  residuals  are  represented  in  these  figures.  The  residuals 

were  computed  with  the  J-B  tables  as  a standard,  and  have  been  averaged  in  blocks  of  5^  in 

azimuth  by  5°  in  angle  of  incidence.  Haw  residuals  were  adjusted  for  station  bias  and  the  Earth's 
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ellipticity  by  using  the  results  in  Sengupta  and  Julian  and  Dziewonski  and  Gilbert,  re$pec- 
tively.  Most  of  these  station  corrections  lie  in  the  range  ±1  sec,  whereas  ellipticity  corrections 
tend  to  make  the  raw  residuals  more  negative  because  the  source  region  is  closer  to  the  equator 
than  to  the  pole. 

Figure  1-3  (a)  shows  a band  of  negative  residuals  that  has  an  orientation  consistent  with  a 
seismic  zone  inclined  between  30°  and  40°  toward  the  west  (the  overall  strike  of  the  seismic 
zone  is  shown  by  the  heavy  bars  in  these  diagrams).  The  minimum  residual  is  —3  sec  arid  the 
average  negative  residual  in  this  band  is  between  —1.5  and  —2.0  sec. 

The  path  length  for  rays  in  the  anomalous  body  can  be  crudely  estimated  from  the  equation 

<aT>  ■ 

where  <AT>  is  the  average  travel-time  residual  corrected  for  all  but  the  source  effect,  fy  is 
the  fractional  velocity  anomaly,  <v>  is  the  average  velocity  for  the  depth  range  d in  the  anom- 
alous body.  Taking  —1.75  sec  to  be  the  corrected  average  residual  from  the  shallow  earthquakes, 
an  average  v^  of  8.5  km/sec  for  the  upper  part  of  the  slab  and  a velocity  anomaly  of  6 percent 
yield  a path  length  of  about  260  km,  which  is  less  than  1/2  the  length  of  the  inclined  seismic 
zone. 

Residuals  from  earthquakes  in  the  depth  range  from  200  to  400  km  [Fig.I-3(b)]  and  from 
depths  greater  than  550  km  [Fig.  1-3 (c))  within  the  same  seismic  zone  show  negative  residuals 
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as  large  as  — 2 sec  that  map  onto  the  same  region  of  the  focal  sphere  occupied  by  the  band  of 
negative  residuals  in  Fig.I-3(a).  The  mean  negative  residual  in  this  region  of  the  focal  sphere 
is  about  —1.0.  An  average  velocity  anomaly  of  +10  percent  in  the  depth  range  550  to  600  km 
can  account  for  a travel-time  residual  of  —1  sec;  however,  this  comparison  of  the  residual  pat- 
tern for  near-surface  and  deep  events  may  not  be  valid  because  of  a steepening  of  the  seismic 
zone  at  depths  greater  than  about  400  km.  Little  evidence  for  such  steepening  can  be  inferred 
from  these  residual  maps.  Estimates  of  absolute  time  residuals  are  required  in  order  to  place 
more  accurate  bounds  on  the  velocity  anomalies  from  deep  earthquakes.  Similar  analyses  will 
be  carried  out  for  the  other  mantle  earthquake  zones  to  show  to  what  extent  the  anomaly  pattern 
. for  the  Tonga-Fiji  region  is  typical  of  such  zones. 

T.  J.  Fitch 
A.  M.  Dziewonski 

D.  IN  SITU  COMPRESSIONAL  VELOCITIES;  REFINEMENTS  TO  THE  INVERSION 
SCHEME  AND  AN  EVALUATION  OF  STATION-SOURCE  BIAS 

Arrival-time  differences  between  a master  and  a number  of  secondary  events  can  be  simul- 
taneously solved  for  the  relative  locations,  in  situ  seismic  velocities,  and  the  direction  cosines 

14 

for  ray  paths  at  the  focus  of  the  master.  In  contrast  to  many  problems  in  geophysics,  this 

one  contains  a set  of  model  parameters  that  is  completely  known  a priori,  although  the  data 

may  be  able  to  resolve  only  a small  subset  of  these  parameters. 

14 

In  Fitch,  direction  cosines  pertaining  to  the  azimuths  of  the  ray  paths  were  taken  to  be 

those  of  the  great-circle  paths,  and  cosines  pertaining  to  the  angles  of  incidence  were  refined 

only  for  paths  that  lay  almost  completely  in  the  upper  mantle.  In  situ  Vp's  5-  to  10-percent 

higher  than  global  average  velocities  for  the  corresponding  depth  ranges  were  computed  for 

three  groups  of  deep  earthquakes.  Ray  paths  that  lay  primarily  in  the  lower  mantle  and  core 

21 

were  defined  by  ray  parameters  computed  from  the  classical  earth  model  of  Jeffreys  or  the 

22 

modern  models  B1  by  Jordan  and  Anderson  and  the  Parameterized  Earth  Model  (PEM)  of 
o 

Dziewonski  et  aL  In  subsequent  work,  the  classical  model  was  used  exclusively  because  the 

bias  imposed  by  a particular  model  was  no  more  than  a few  percent  in  velocity. 

1 5 

Fitch  showed  that  in  the  depth  range  from  450  to  650  in  the  Tonga  seismic  zone,  m situ 

V is  anomalously  high  only  at  depths  greater  than  550  to  600  km.  These  results  have  now  been 

P 1 0 23  24 

confirmed  with  improved  data  sets  and  damped  least-square  solutions  ’ ’ akin  to  a sto- 

25  26 

chastic  inverse  proposed  by  Franklin  and  Jordan  and  Franklin. 

The  P arrival  times  were  taken  from  the  Bulletin  of  the  ISC.  Chi-square  tests  of  the  re- 
siduals of  the  differential  times  are  consistent  with  standard  deviations  of  about  0.6  sec,  when 
grossly  incorrect  times  are  eliminated.  The  corresponding  errors  in  the  individual  arrival 
times  are  about  0.4  sec.  Solutions  with  synthetic  data  consistent  with  the  J-B  travel  times, 
actual  station-source  configurations,  and  the  above-mentioned  level  of  uncertainty  show  that 
the  separations  between  master  and  secondary  events  should  be  no  less  than  about  30  km  and 
no  more  than  about  100  km  to  be  reasonably  confident  that  the  computed  in  situ  V^  lies  within 
two  standard  deviations  of  the  correct  value.  The  latest  compilation  of  P times  from  the  Tonga 
seismic  zone  (carried  out  while  the  author  was  still  at  the  Australian  National  University)  re- 
duced the  largest  separation  between  master  and  secondary  earthquakes  to  about  100  km  for 
each  of  the  stable  solutions  in  Table  1-1. 
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TABLE  1-1 


COMPARISON  OF  OBSERVED  AND  CALCULATED  VELOCITIES 


Depth  of 

Master  Earthquake 
(km) 

Largest 
Separation 
from  Master 
(km) 

Vp  from 
Synthetic  Data 
(km/sec) 

Vp  from 
Real  Data 
(km/sec) 

Vp  from 
Jeffreys^ 
(km/sec) 

493 

102 

9. 4 ± 0.  2 

9.  7 ± 0.  2 • 

9.60 

520 

92 

Unstable 

Unstable 

9.76 

535 

75 

10.1  ±0.2 

9.4  ± 0.*2 
(marginally  stable) 

9.88 

572 

126 

9.  8 ± 0. 2 

Unstable 

10.09 

597 

93 

9.9  ±0.2 

10.7  ±0.2 

10.  23 

608 

144 

Unstable 

Unstable 

10. 29 

624 

107 

10.2  ±0.2 

10.7  ±0.2 

10.37 

650 

91 

9.  8 ± 0.  2 

11.2  ± 0.2 

10. 47 

These  solutions  include  refined  estimates  of  the  azimuths  as  well  as  the  angles  of  incidence 

for  ray  paths  to  the  stations  shown  in  Fig.  1-4.  Rather  than  a simple  least-squares  solution  that 

minimizes  |b  — Ay|2  where  x and  b are  the  data  and  model  vectors,  respectively,  and  A is  the 

matrix  of  partial  derivatives,  these  solutions  were  damped  by  the  addition  of  a small  component 

12  2 2 

to  the  diagonal  elements  of  the  solution  matrix.  This  is  equivalent  to  minimizing  |b  — Axj  + 0 x . 
where  02  is  the  vector  containing  the  damping  coefficients.  Following  Aki  et  al.10  thesp  coeffi- 


2 2 ^ 
cients  were  initially  set  equal  to  v/v  , where  vc 

7 * 

data, 


is  the  variance  of  the  random  errors  in  the 


and  the  a are  estimates  of  the  variance  in  the  model  parameters.  To  achieve  mbre 
rapid  convergence  to  a solution,  the  initial  damping  coefficients  were  multiplied  by  0.1. 

Damped  solutions  are  preferred  to  the  simple  solutions  because  more  rigorous  account  is 
taken  of  errors  in  the  data  that  inevitably  limit  the  resolution  of  the  model  parameters.  The 
least-resolved  parameters  are  the  angles  defining  the  direction  cosines  and  the  separation  be- 
tween events.  The  in  situ  velocity  in  all  cases  was  resolved  to  greater  than  95  percent* 

Solutions  with  synthetic  data  reveal  a bias  in  the  in  situ  velocities  that  is  imposed  by  the 

particular  station-source  configuration  in  the  real  earth.  This  bias  was  thought  to  be  negligible 
15 

in  Fitch;  however,  the  results  in  Table  1-1  and  Fig.  1-5  show  that  this  bias  could  be  a3  great 
as  5 percent  for  these  particular  data  sets.  The  in  situ  velocities  pertaining  to  the  master 
earthquakes  at  depths  of  650  and  597  km  are  significantly  less  than  the  corresponding  Jeffreys; 
velocities  in  the  solutions  with  synthetic  data,  and  greater  in  the  solutions  with  real  data.  The 
bias  in  the  solutions  with  synthetic  data  results  from  the  locations  of  these  master  earthquakes 
which  are  below  all  of  their  secondary  events  (Fig.  1-5). 

The  higher  velocities  computed  with  the  real  data  cannot  be  explained  by  station-spurce  bias. 

14  15 

Consequently  these  solutions,  as  well  as  less-sophisticated  solutions  reported  by  Fitch|,  ’ 
reveal  a positive  velocity  anomaly  near  the  base  of  this  mantle  seismic  zone.  The  highest 
in  situ  velocity  was  computed  for  a group  of  events  with  the  master  at  650  km,  which  iP 
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TABLE  1-2 

STRUCTURE-DEPENDENT  RAY  PATHSt 

Azimuth  (deg) 

Angles  of  Incidence  (deg) 

JB 

Refined 

Jeffreys21 

Refined 

AFI 

493 

37.2 

39.2  ± 4.7  (0.859) 

75.2 

75. 3 ± 5. 3 (0.  765) 

597 

58.3 

58.8  ± 6.0  ( 0.830) 

109.6 

128.2  ± 5.8  (0.903) 

624 

60.3 

43. 1 ± 3. 1 (0. 940) 

107.7 

121. 1 ± 5.4  (0.618) 

650 

46. 1 

46.0  ± 4.9  (0.873) 

96.8  * 

109.5  ± 4.6  (0.876) 

GNZ 

493 

186.4 

189.4  ±6.6  (0.710) 

73.8 

79. 3 ± 5.  7 (0.  705) 

597 

187.9 

183.6  ± 6.7  (0.552) 

67.4 

66. 8 ± 6. 5 (0.  709) 

624 

186.4 

199.  1 ± 3.3  (0.767) 

69.2 

66.2  ± 7. 1 (0.907) 

WEL 

493 

193.5 

202.5  ±6.5  (0.981) 

66.9 

69.  7 ± 5. 4 (0. 952) 

624 

191.5 

197. 1 ± 3.  7 (0.  779) 

65.8 

60. 9 ± 7. 5 (0. 924) 

650 

192.9 

187.8  ± 5. 1 (0.236) 

64.0 

83.2  ± 5.0  ( 0.773) 

RAO 

493 

157.5 

156.0  ± 4.8  (0.886) 

122.6 

109.2  ±5.8  (0.851) 

597 

177.8 

162.0  ± 7.9  (0.894) 

92.7 

88.3  ± 5.8  (0.947) 

650 

172.3 

185.0  ± 5.8  (0.876) 

112.2 

124.6  ± 5.2  (0.955) 

t Corresponding  diagonal  element  of  the  resolution  matrix  is  given  in  parentheses. 
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nominally  the  depth  of  a major  transition  zone  between  the  upper  and  lower  mantles.  The  PEM 
velocities  in  Fig.  1-5  illustrate  the  size  of  the  velocity  contrast  that  is  consistent  with  a combina- 
tion of  seismological  and  petrological  prejudice. 

By  definition,  deep  earthquakes  occur  in  anomalous  mantle.  A map  view  of  the  Tongb  deep 
seismic  zone  is  shown  in  Fig.  1-6  along  with  the  epicenters  of  the  master  earthquakes.  A signif- 
icant fraction  of  the  paths  to  the  stations  AFI,  RAO,  GNZ,  and  WEL  (Fig.  1-4)  lie  within  dr 
closely  parallel  the  seismic  zone  and  thus  are  expected  to  show  the  greatest  deviations  from 
the  paths  predicted  by  symmetrical  earth  models. 

In  Table  1-2  the  refined  azimuths  that  differ  from  the  predicted  azimuths  by  more  th^n  one 
standard  deviation  (underscored  in  the  table)  reveal  rotations  of  ray  paths  toward  the  stride  of 
the  seismic  zone.  (The  azimuth  for  the  path  WEL650  is  an  apparent  exception;  however,  it  is 
essentially  unresolved  by  the  data.)  Similarly  the  angles  of  incidence  from  the  deeper  master 
earthquakes  to  AFI  (Fig.  1-4)  are  rotated  toward  the  upward  vertical,  and  thus  a larger  frjaction 
of  the  corresponding  ray  paths  will  lie  in  the  seismic  zone. 

Taken  together,  the  in  situ  velocities  and  the  direction  cosines  for  ray  paths  to  the  closer 
stations  are  consistent  with  other  seismological  and  geophysical  evidence  for  a body  of  anomalous 
material  containing  mantle  earthquakes.  However,  the  in  situ  velocities  show  that  the  vqlocity 
anomaly  associated  with  this  body  is  not  uniformly  distributed  in  the  depth  range  from  450  to 

600  km.  One  way  to  explain  the  anomaly  pattern  is  by  an  elevated  650-km  discontinuity  qs 
1 5 

proposed  in  Fitch.  ^ j 

E.  IN  SITU  RATIOS  OF  COMPRESSIONAL-TO-SHEAR  VELOCITY 
IN  SHALLOW-  AND  DEEP-SOURCE  REGIONS 

The  seismic  velocities  V and  V and  their  ratio  are  being  reported  in  many  studies  both 

P s 

within  and  outside  the  Earth's  major  earthquake  zones.  The  inferred  velocities  or  velocity 

ratios  are,  in  almost  every  case,  based  on  less  than  ten  arrival  times  of  the  P and  S phases 

from  a single  event.  These  inferences  from  such  a limited  data  set  are  vulnerable  to  various 

pathologies  that  may  exist  within  an  active  fault  zone  or  between  the  source  region  and  the  re- 
28 

ceivers.  For  example,  anomalously  high  attenuation  of  weak  signals  may  lead  to  arrival 
times  which  are  picked  systematically  late,  and  thus  anomalously  low  velocities,  or  S-tq-P 
conversions,  may  be  misinterpreted  as  the  S phase,  in  which  case  the  inferred  shear  velocity 
will  be  too  fast. 

29 

A major  advantage  of  the  inversion  scheme  proposed  by  Fitch  and  Rynn  and  promoted 

here  is  that  the  velocity  ratio  is  computed  from  a large  data  set.  Several  hundred  arrival 

times  from  10  to  15  events  are  simultaneously  inverted  for  relative  locations  and  elevation 

angles  for  ray  paths,  as  well  as  the  velocity  ratio.  With  such  a data  set,  large  reading  (errors 

contained  in  a subset  of  the  data  are  likely  to  be  discovered.  Other  advantages  of  this  scheme 

are  that  the  velocity  ratio  pertains  to  the  source  region  and,  in  principle,  does  not  require  any 

knowledge  of  the  velocity  structure  outside  that  region. 

29 

The  inversion  scheme  reported  by  Fitch  and  Rynn  is  generalized  so  that  P and  S waves 
traveling  between  each  source  and  receiver  are  represented  by  independent  ray  paths.  Resolu- 


tion must  be  sacrificed  to  achieve  this  generalization  for  the  case  of  data  from  crustal  earth- 
quakes that  are  only  locally  recorded.  The  loss  of  resolution  occurs  mainly  in  the  rel 
depths  and  in  the  direction  cosines  for  ray  paths  at  the  focus  of  the  master  earthquake 


V /v 

p'  s 
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TABLE  1-3 

A SAMPLE  SOLUTION 


Secondary 

Event 

No. 

Differential 
Origin  Time 
(sec) 

Relative  Locations  of  Secondary  Events 

Distance 

(km) 

Azimuth 

(deg) 

Elevation  Angle 
(deg) 

1 

10.04  ± 0.03 

1.06  ± 0.  15 

41.64  ± 4.58 

62. 62  ± 8.  59 

2 

10.02  ±0.03 

0.87  ± 0.  16 

27.35  ± 5.73 

58.  74  ± 8.  02 

3 

10.47  ± 0.02 

1.83  ± 0.  14 

174.81  ± 4.58 

142.33  ± 2.86 

4 

10. 15  ± 0.02 

1.05  ± 0.09 

344.86  ± 6.30 

108.42  ± 8.02 

5 

10.49  ± 0.02 

1.77  ± 0.  14 

171.79  ± 4.01 

133.92  ± 2.86 

6 

10.02  ± 0.02 

2.04  ± 0.  15 

30.22  ±3.44 

53.69  ±3.44 

7 

10. 11  ±0.03 

1.47  ± 0. 12 

2.51  ± 4.01 

83.56  ± 6.30 

8 

10.33  ± 0.02 

1.06  ± 0.  13 

146.50  ± 7.45 

135.95  ± 4.58 

9 

10.02  ± 0.02 

1.65  ± 0.  14 

30.39  ± 2.86 

58.28  ±3.44 

10 

9.99  ± 0.02 

1.06  ± 0.  11 

116.50  ± 6.88 

36. 18  ± 5. 16 

11 

10.27  ± 0.02 

1.08  ± 0.07 

39.  15  ± 5.  16 

131. 16  ± 8.59 

12 

10.45  ± 0.02 

1.55  ± 0.  11 

166.02  ± 4.01 

117.93  ± 4.01 

13 

9.99  ± 0.02 

2.02  ± 0. 14 

25.  78  ± 2.  86 

52.  88  ± 2.  86 

Station 

Ray-Path  Elevation  Angles^ 

Initial  VD/V$ 
(deg) 

Refined  V (deg) 

P Wave 

S Wave 

LIV 

168 

165.88  ± 3.95 

170.54  ± 4.52 

TAB 

156 

155.84  ± 2.87 

166.16  ±3.04 

DAM 

155 

157.81  ± 2.67 

158.64  ± 2.80 

ORV 

154 

155.05  ± 2.02 

156.45  ±2.18 

BID 

143 

141.86  ± 5.00 

- 

WYN 

140 

141.61  ± 7.34 

122.23  ± 6.57 

LON 

136 

120.38  ± 5.07 

- 

CAM 

137 

133.55  ± 4.46 

- 

FOR 

131 

141.  10  ± 2.03 

- 

RAT 

130 

136.74  ± 4.52 

126.29  ± 4.32 

FIG 

127 

128.  15  ± 5.45 

124. 12  ± 8. 19 

KAT 

120 

128.92  ± 4.46 

- 

HON 

114 

106.07  ± 6.22 

- 

STI 

118 

1 04.  10  ± 4.36 

- 

Initial  y/2. 

Refined  1.60  ± 0.  02  with  the  corresponding  element  of  the  resolution  matrix  0.999. 
X2=  105  and  <x2  > = 149. 


71  of  74  eigenvalues  retained. 
. t 0°  is  the  downward  vertical. 
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is  nearly  perfectly  resolved  in  solutions  with  real  and  synthetic  data  of  high  quality,  i.e.,,  arrival 
times  with  random  errors  no  larger  than,  say,  0.05  sec.  Several  data  sets,  each  containing 
between  200  and  300  arrival  times,  were  compiled  from  the  better-recorded  earthquakes  in  the 
NW  corner  of  the  aftershock  zone  of  the  1 August  1975  Oroville,  California  earthquake. 

The  preferred  estimate  of  in  situ  V /V  lies  in  the  range  1.65  to  1.60,  which  is  significantly 

P s 

less  than  1.70  — the  ratio  assumed  by  the  Office  of  Earthquake  Research  at  Menlo  Park,  California 
for  computing  the  absolute  locations  of  these  events.  The  solution  summarized  in  Ts.b\e  1-3  is 
typical.  The  data  for  this  solution  come  from  a source  region  that  is  annular  in  shape,  with  the 
master  located  in  the  center.  The  inner  radius  of  about  1 km  is  limited* by  the  errors  ifi  reading 
arrival  times,  and  the  outer  radius  is  limited  by  the  condition  that  the  separation  between  master 
and  secondary  events  be  small  by  comparison  with  the  length  of  the  ray  path  to  the  nearest  sta- 
tion (Fig.  1-7).  The  standard  errors  in  high-quality  P and  S times  were  found  to  be  between 
±0.01  to  ±0.02  sec  and  between  ±0.03  to  ±0.04  sec,  respectively.  These  levels  of  uncertainty 
are  consistent  with  chi-square  tests  on  residuals  of  the  differential  P and  S times. 

The  relative  locations  projected  on  a focal  mechanism  solution  for  the  main  shock  (pig.  1-8) 
define  a broad  band  from  north  to  south  between  the  nodal  planes.  The  relative  locations  pertain 
to  activity  near  the  base  of  the  aftershock  zone  at  a depth  between  8 and  10  km.  Apparently, 
in  this  depth  range  the  seismic  zone  has  nearly  the  same  strike  as  shallower  parts  of  the  zone, 

but  is  inclined  nearly  vertically.  This  can  also  be  seen  in  the  vertical  distribution  of  aijter- 

30 

shocks  reported  by  Lahr  et  al. 

The  stability  of  the  inversion  scheme  increases  as  the  distance  between  master  and  secondary 
events  increases.  However,  these  distances  must  be  sufficiently  small  to  avoid  making  assump- 
tions about  the  geometry  of  the  ray  paths  outside  the  source  region.  Larger  source  regions  are 
permitted  at  greater  depth;  consequently,  the  inversion  scheme  becomes  a more  robust  esti- 
mator of  the  velocity  ratio  as  source  depth  increases.  This  is  demonstrated  by  solutions  with 
synthetic  data  from  a hypothetical  source  region  centered  at  a depth  of  100  km. 

Computations  reported  here  were  carried  out  at  the  computer  center  of  the  Australian 
National  University.  T j pitch 

A.  G.  Lindh^ 

F.  A MOMENT-TENSOR  REPRESENTATION  OF  BODY-WAVE  DISPLACEMENT 

VECTORS  ON  THE  FOCAL  SPHERE 

Presently  two  point-source  moment -tensor  representation  theorems  exist;  the  original 
free  oscillation  expansion  by  Gilbert,^1  and  the  Love-  and  Rayleigh- wave  expansion  done  by 
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McCowan/4-  The  advantages  to  the  moment-tensor  approach  are  well  known.  These  ex|>ansions 
represent  nonpropagating  point  sources  in  their  most  general  form,  including  all  three  {possible 
mechanisms;  explosive,  double  couple,  and  compensated  linear  vector  dipole.  Furthermore, 
the  resulting  vector  displacement  is  always  a linear  combination  of  the  moment-tensor  t demerits 
thereby  making  the  determination  of  an  unknown  source  mechanism  amenable  to  linear -estimation 
theory.  The  purpose  here  is  to  extend  the  moment-tensor  formulation  to  body-wave  displace- 
ment vectors  on  the  focal  sphere. 

The  problem  is  solved  in  spherical  coordinates  using  the  geometry  shown  in  Fig.  1-9.  The 
source  is  assumed  to  act  at  the  center  of  a homogeneous  elastic  solid  where  it  is  described  by 

A A * 

the  fixed  Cartesian  coordinate  system  (x,y,z)  with  associated  unit  vectors  (i,  j,k).  The  (Jisplace- 
ment  vector  resulting  from  this  source  at  a distance  r from  the  origin  oriented  at  the  a,ngles 


t National  Center  for  Earthquake  Research,  USGS,  Menlo  Park,  CA. 
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0 and  <p  from  the  source  coordinate  system  will  be  represented  by  its  radial,  polar,  and 

azimuthal  components:  s . s^,  and  s M with  respect  to  the  unit  vectors  £ , £^,  £ . Then  the 

r 0 <p  r Q (p 

usual  arguments  lead  to  the  point-source  moment -tensor  solution: 


1 

0 

0 


s*(r) 


f 9G(r,r ')  8G(r,r')  8G(r\P)l 
l 8x»  dyr  3z'  Jrt  = 0 


0 

1 

0 

0 


0 

1 


d-1) 


Here,  s (r  ) is  the  vector  of  the  three  displacement  components  mentioned  above,  G(r  ,r')  is  the 


Green's  tensor  in  £ , £_ 
r 0 

k components 


tensor  in  i , j , 


£ ^ components  (Morse  and  Feshbach 


33 


),  M is  the  source  moment 


and  A is  the  following  tensor: 

'sin0  cos  <p  sin0  sirup 

cos  0 cos  <p  cos0  sirup 

— sirup  cos  (p 


d-2) 


The  last  column  vector  is  necessary  to  perform  the  remaining  indicial  sum. 

33 

As  shown  by  Morse  and  Feshbach  (Chap.  13),  the  Green's  tensor  can  be  decomposed  into 


longitudinal  and  transverse  parts  which  propagate  P and  S waves,  respectively, 
infinite  domain  expressions  for  these  are,  respectively, 


The  correct 


1 f/1+ika 

pa2  [\  k2R2 

■*K- 


R 


I - 


RRT 

~rT 


+ 3ik  R-k2R2Y| 

a a \\ 

k 2R2  / 


-ik  R 

% ot 


ikpR  -kg2R2 


--T  / 3 + 3ik.  R — k2R2 


I + 


RR 

R2 


P 


)] 


;V 

R 


(1-3) 


It  is  easiest  to  refer  to  Morse  and  Feshbach  for  an  explanation  of  the  symbols  and  conventions 
in  these  expressions.  The  solutions  can  now  be  obtained  by  substituting  L and  T from 
Eqs.  (1-3)  for  G in  Eq.  (I- 1 ) and  performing  the  indicated  operations. 

After  discarding  the  terms  which  decay  faster  than  l/r,  the  following  are  obtained: 


sp(r,t) 

sg(r,t) 


P(M,0,q 


P f+0°icpg(WRiw(t-r/o)dwfr 


2i rp  or  * -oo 


U—  C icog(w)fiw(t'r/'3)dW[S0(M.0,<?>)  i + S (M.0.*)  i ] 


d-4) 
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Here  a represents  the  P-wave  velocity,  p the  shear-wave  velocity,  p the  density,  and  g(a;)  the 
source- excitation  function.  The  angular  functions  are: 

2 2 2 2 2 
P(M,0  ,<p)  = M sin  0 cos  <p  + M sin  0 sin  q>  + M cos  0 
^ xx  yy  zz 

+ M sin20  sin2<p  + M sin  20  cos  w + M sin  20  sin<p 
xy  Y xz  T yz  Y 

M ? M ? M 

Sq{M, Q,<p)  - -£X-  sin  20  cos  <p  + — ^ sin  20  sin  <p — ^ sin  20 


M 

+ 2^  sin  20  sin  Zip  + M cos  20  cos  <p  + M cos  20  sin  q 

c xy  yz 


M M 

S (M ,0,<p)  =-^*sin0  sin 2 w + — ^sin0  sin2<p  + M *0 

<p  Y 2 T 2 ^ zz 

+ M sin0  cos  2 tp  — M cos  0 sin<p  + M cos  0 cos  <p 
xy  xz  yz 

For  the  case  of  a step-source  time  function  g(c o)  = — l/icj,  Eqs.  (1-4)  reduce  to: 

Sp(r.t)  = 6 (t  - r/a)  f r 

pa  r 

s"o(r,t)  = 6(t  ~('/g)  [S  (M,e,<p)  / + S (M,e,<p)  i ) . 


(1-5) 


(1-6) 


As  can  be  seen  in  Eqs.  (1-5),  an  isotropic  moment  tensor  M = ml  produces  monopqle  radia- 
tion. Furthermore,  a double-couple  moment  tensor 


M = m(pd  T + dp 


(1-7) 


where  p and  d are  the  normal  and  slip  unit  vectors,  respectively,  yields  solutions  wllich  agree 

34 

with  those  of  Ben-Menahem  and  Singh.  Equations  (1-5),  however,  are  the  most  general  solu- 
tion to  the  point-source  problem  because  they  include  all  three  possible  source  mechanisms. 
Extended  sources  may  be  constructed  by  integrating  these  solutions  over  the  source  volume. 

D.  W.  McCowan 
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Fig.I-1.  Short-period  waveforms  recorded  at  SRO  stations.  Each  seismogram  is  identified 
by  station  (A  = AMNO,  M = MAIO,  G = GUMO),  depth  h (km),  and  distance  A (deg).  Maxi- 
mum amplitude  of  signal  is  also  given.  Each  seismogram  is  of  20-sec  duration,  centered  at 
arrival  time  predicted  by  J-B  tables. 


89  m/i 


38  m/i 


194  m/i 


477  m/i 


(a) 


Fig.  I-2(a).  Azimuthal  dependence  of  P travel-time  anomalies  for  western  and  central  U.S., 
as  referred  to  model  PEM,  for  distance  ranges  A = 30°  to  40°.  Azimuth  range  is  divided 
nirtor 45 * regions : —flanks indicate  no;  or  not"enongh7~rlata; — * 


LATITUDE  <#N) 


(b) 


Fig.  1-2 (b).  Azimuthal  dependence  of  P travel-time  anomalies  for  western  and  central  U.S., 
as  referred  to  model  PEM,  for  distance  ranges  A = 80°  to  90°.  Azimuth  range  is  divided 
into  45°  regions.  Blanks  indicate  no,  or  not  enough,  data. 


Fig.  1-3.  Residual  spheres.  Residuals  in  seconds  are  mapped  in  an  equal  area  sense 
onto  lower  half  of  focal  sphere.  Data  have  been  taken  from  events  ,with  latitude  riange 
15°S  to  35 °S  and  longitudes  greater  than  177  °E  in  Tonga-Kermadec  region,  and  station 
corrections  due  to  Sengupta  and  Julian* 9 have  been  included.  Strike  of  seismic  zone  is 
shown  by  heavy  line  oriented  N23°E.  Depth  ranges  shown  are  (a)  depth  <100  km, 
(b)  200  < depth  < 400  km,  and  (c)  depth  > 550  km. 
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Fig.  1-4.  Southwest  Pacific  and  ray  paths  to  island  stations. 


8 9 10  it 

Vp  (km/sec) 


Fig.  1-5.  Vp  vs  depth.  Error  bars  on  in  situ  V^’s  are  one  standard  deviation. 

Vertical  line  attached  to  each  synthetic  datum  gives  depth  range  for  particular 
group  of  secondary  events. 
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Fig.  1-6.  Deep  seismic  zone:  a map 

view.  Dashed  line  indicates  location 
of  epicenters  of  earthquakes  with  depth 
500  km  (after  Sykes27). 
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Fig.  1-7.  Station-source  configuration 
for  real  and  synthetic  data. 
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Fig.  1-8.  Relative  locations  projected  on  focal  mechanism  for  main  shock. 
Nodal  planes  and  relative  locations  are  projected  in  an  equal  area  sense 
onto  lower  half  of  focal  sphere.  L.  refers  to  separation  between  master 
and  jth  secondary  events.  ^ 
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II.  SURFACE-WAVE  STUDIES 


A.  AN  EFFICIENT  ALGORITHM  FOR  DETERMINING  PHASE  AND  GROUP 
VELOCITIES  OF  SURFACE  SEISMIC  WAVES 

The  method  of  bandpass  filtering  in  the  frequency  domain  has  become  a powerful  topi  for 

the  determination  of  phase-  and  group-velocity  plots  of  surface  seismic  waves.  The  original 

seismic  signal  is  first  filtered  to  remove  all  frequencies  except  a narrow  positive  frequency 

band  centered  about  some  frequency  co The  group  arrival  time  at  frequency  cx)Q  is  then  defined 

as  the  time  of  maximum  amplitude  in  the  envelope  of  the  complex  filtered  waveform,  and  the 

phase  velocity  can  be  determined  from  the  phase  relationship  at  the  time  of  maximum  amplitude. 

The  process  is  repeated  for  several  center  frequencies  spanning  the  range  from  approximately 

0.007  to  0.06  Hz  to  construct  group-  and  phase -velocity  curves.  The  method  is  particularly 

useful  in  situations  when  the  signal  is  not  well-dispersed,  as  in  typical  oceanic  paths  or  short 

1-4 

continental  paths,  for  in  these  cases  a simple  peak  and  trough  method  is  impossible. 

Since  we  have  available  to  us  an  extensive  database  of  high-quality  digital  SRO  surface  - 
wave  data,  it  would  be  worthwhile  to  construct  group-  and  phase-velocity  curves  for  a large 
number  of  paths,  both  continental  and  oceanic.  The  limitation  of  the  method  is  that  it  requires 
an  excessive  amount  of  computation  time,  and  since  we  plan  to  run  a large  data  set  through  the 
processing,  it  is  advantageous  to  reduce  the  computation  time,  if  at  all  possible.  Therefore, 
we  have  developed  a slight  modification  to  the  original  algorithm  which  results  in  a greater 
than  10:1  savings  in  computation  time. 

In  the  original  algorithm,  the  Fourier  transform  of  the  original  waveform  is  multiplied  by 

the  Fourier  transform  of  the  bandpass  filter,  and  an  inverse  Fourier  transform  of  the  product 

yields  the  filtered  seismic  signal.  The  FFT  size  must  be  large  to  obtain  adequate  frequency 

resolution.  (We  used  a 2048-point  FFT  for  data  sampled  at  2-sec  intervals  to  obtain  0.00025-Hz 

resolution. ) The  spectrum  of  the  filter  is  zero  everywhere  except  in  a narrow  band  about  go  . 

1 c 

We  used  the  Gaussian  filter  suggested  by  Dziewonski: 


H(w)  = exp  {-a  [(a> 

where  a is  typically  50.  However,  our  proposed  modification  is  suitable  for  any  narr 
pass  filter. 

Since  the  spectrum  of  the  filtered  seismogram  is  zero  outside  of  a narrow  frequency  range, 
a considerable  reduction  in  bandwidth  can  be  achieved  by  simply  shifting  the  spectrum  by  an 
amount  —o>c  (Fig.  II- 1 ).  A frequency  shift  of  — o>c  is  equivalent  to  multiplication  of  each  sample 
s(n)  of  the  time-domain  signal  by  exp  [— jco  nt],  where  t is  the  sampling  interval.  The  advantage 
of  such  a step  is  that  the  size  of  the  IFFT  can  now  be  reduced  to  include  only  the  nonzero  sam- 
ples near  c without  effecting  any  loss  of  information.  Whereas  a 2048 -point  IFFT  previously 
was  needed  for  each  u >c,  now  only  a 128-  or,  at  most,  256-point  IFFT  is  needed.  The  large 
forward  FFT  is  computed  only  once,  whereas  the  IFFT  must  be  computed  once  for  each  center 
frequency.  We  have  found  it  necessary  to  include  75  center  frequencies  'in  order  to  resolve 
phase  ambiguities;  hence,  the  forward  FFT  time  is  insignificant  by  comparison. 

The  output  of  the  modified  Dziewonski  algorithm  is  a signal  with  slowly  varying  amjplitude 
and  slowly  varying  phase  relationship  (Fig.  II-2).  Group  velocity  is  computed  as: 
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Ug  =r/(delay  + npkt-t.nst) 


where  r is  the  great-circle  distance  between  earthquake  epicenter  and  receiver,  delay  is  the 
time  between  the  event  and  the  first  sample  of  the  seismogram,  n^t  is  the  time  between  the 
origin  and  the  peak  amplitude  in  the  filtered  seismogram,  and  tngt  is  the  instrument  time  delay. 
The  peak  amplitude  location  is  determined  by  matching  a parabola  to  the  three  samples  surround- 
ing the  peak  and  choosing  the  interpolated  peak  time  as  follows: 

y^  = A(n  — 1)  — A(n)  n = peak  location 

y3  = A(n  + 1)  — A(n)  A(n)  = amplitude  at  peak 

(yt  -y3) 

pk  2(yi+y3) 

Figure  II-3  shows  a typical  group-velocity  curve  for  a well-dispersed  continental  path. 

The  computation  of  phase  velocity  is  slightly  more  complicated.  The  phase  at  peak  ampli- 
tude in  the  original  Dziewonski  algorithm  is  defined  approximately  as  (ignoring,  for  the  moment, 
instrument  and  source  phase): 

«’c<npk)  = “cV*  + Wc  * delay  - kcr  * 2?ri 

where  k^  is  the  wavenumber  corresponding  to  frequency  cjc,  and  Zni  is  an  unknown  number  of 
full  cycles.  In  the  modified  algorithm,  the  first  term,  Gt^n^t,  is  removed  by  the  frequency 
shift  of  wc,  and  hence  the  equation  is  somewhat  simplified.  The  first  step  is  to  compile  a se- 
quence of  phase  values  corresppnding  to  maximum  envelope  amplitude  for  75  consecutive  fre- 
quencies o^.  Then,  this  sequence  is  unraveled  to  resolve  ambiguities  of  27T  (made  continuous) 
and  corrected  for  instrument  response,  and,  if  possible,  source  phase.  Next,  cfa;^)  is  evalu- 
ated at  a low  frequency  for  several  values  of  i: 

<Pc(npk)  = wc  ' delay  ~ kcr  * 2ni 

& u r w r 

c c c 

C(Wc>  = \ - V = wc  ' delay  ~ <z>c(npk)  ± 2?ri 

The  first  reasonable  choice  for  i (e.g.,  velocity  less  than  4.6)  is  kept,  and  phase  velocities  are 
then  computed  for  all  75  frequencies  using  i,  i + 1,  and  i + 2 to  construct  three  possible  phase - 
velocity  curves.  The  correct  curve  can  then  usually  be  selected  by  eye,  given  a knowledge  of 
what  is  reasonable  (Fig.  II - 3 ) . g Seneff 

B.  GROUP- VELOCITY  MEASUREMENTS  ACROSS  EURASIA  FROM  MASHAD  SRO 

During  January- March  1976,  some  100  events  within  or  on  the  margins  of  Eurasia  were 
sufficiently  well-recorded  at  Mashad  SRO  such  that  well-dispersed  surface-wave  trains  were 
observable.  The  wide  dynamic  range  of  the  instrument  is  particularly  well-suited  to  dispersion 
studies  since  small  events  from  regions  of  low  seismicity  produce  usable  data.  The  three- 
component  data  have  been  rotated  into  the  source  azimuth  and  this,  in  general,  produces  good 
separation  of  Love-  and  Rayleigh-wave  radiation  into  the  transverse  and  radial  components. 

To  date,  only  the  vertical  (Rayleigh)  component  has  been  studied;  in  the  future,  we  intend  to 
utilize  the  polarization  properties  of  surface  waves  to  reduce  body-wave  contamination. 
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Figure  II-4  shows  the  events  to  which  Rayleigh-wave  dispersion  has  been  calculated  from 
Mashad.  Many  more  events  during  the  remainder  of  1976,  giving  more  complete  coverage  in 
distance  and  azimuth,  remain  to  be  analyzed.  A new  adaption  of  the  multiple -filtering  tech- 
nique^ (see  Sec.  A above),  permitting  rapid  and  simultaneous  determination  of  group  and  phase 
velocities,  has  been  used  to  produce  dispersion  curves  for  these  paths. 

Although  the  phase-velocity  resolution  is  much  better  than  that  in  group  velocity,  the  small- 
ness of  many  of  the  events  means  that  the  crucial  initial  source  phase  is  virtually  impossible  to 
determine  because  of  the  difficulty  of  obtaining  both  source  mechanism  and  depth.  For  fevents 
in  regions  where  the  relative  motion  of  the  plates,  and  thus  the  dominant  slip  vector,  is  known 
(such  as  the  western  boundary  of  the  Pacific  plate),  it  is  possible  to  make  an  intelligent  guess 
as  to  the  source  mechanism;  but,  within  continental  Asia,  the  variety  of  observed  faultihg  even 
within  small  source  regions  makes  this  very  difficult. 

Figure  11-5  shows  group-velocity  dispersion  determined  for  some  of  the  paths  showh  in 
Fig.  II -4.  By  far,  the  most  notable  effect  is  that  of  the  very  thick  low- velocity  crust  in  Tibet 
and  Western  China,  which  lowers  group  velocities  to  below  3.0  km/sec  even  at  periods  of  40  to 
50  sec.  Dispersion  to  Siberia,  Mongolia,  and  the  Kuriles  is  somewhat  faster  than  expected, 
considering  the  varied  tectonic  regimes  crossed.  The  highest  group  velocities  are  those  to 
Iceland  and  Svalbard  in  the  north-west  quadrant,  crossing  the  Russian  platform  and  Baltic 
Shield. 

Many  more  events  remain  to  be  analyzed,  and  when  sufficient  data  have  been  obtained  we 
intend  to  carry  out  an  inversion  to  determine  group  velocity  as  a function  of  geographical  posi- 
tion within  the  region,  without  making  any  prior  assumptions  as  to  the  connection  between  sur- 
face geology  and  dispersion.  R North 


C.  GROUP-VELOCITY  DISPERSION  MEASUREMENTS  AT  NEAR-IN  DISTANCES 

Of  the  three  commonly  used  methods  for  measuring  the  group-velocity  dispersion  of  a 
wavetrain  — instantaneous  frequency  analysis,  narrowband  filtering,  and  moving-window  spec- 
tral computation  — only  the  latter  is  easily  adaptable  for  use  on  compact  waveforms.  The  in- 
herent weakness  in  the  moving -window  method  has  always  been  obtaining  adequate  spectral 
bandwidth  from  short  segments  of  data.  If  the  data  sample  is  long,  the  computed  spectrum  is 
characteristic  of  the  whole  window,  and,  if  the  data  sample  is  short,  it  is  difficult  to  resolve 

low  frequencies.  However,  recently,  the  addition  of  Maximum  Entropy  Spectral  Analysis 
6 7 

(MESA)  to  data-adaptive  prediction-error  filtering  offers  a potentially  valuable  means  of 
achieving  adequate  frequency  resolution  from  short  data  windows. 

The  method  consists  of  computing  the  data-adaptive  prediction-error  operator  using 
7 8 

Griffiths  formulation  of  the  Widrow  algorithm.  This  is  basically  the  same  algorithm  used 
for  SP  deconvolution  studies  reported  in  the  previous  two  SATS.  At  specified  group  arrival 
times,  the  inverse  of  the  prediction-error-operator  spectrum  is  computed  and  contour-plotted. 
The  advantage  to  using  an  adaptive  filtering  algorithm  is  that  the  data  design  the  filter,  Ithus 
allowing  the  filter  to  “track"  variations  in  the  data  spectrum. 

As  examples  of  the  method,  group-velocity-dispersion  curves  were  computed  for  two  rela- 
tively close-in  events  selected  from  our  SRO  data  catalog.  These  data  are  particularly  well- 
suited  to  dispersion  computations  because  of  their  broad-band  high- dynamic -range  chartacter- 
istics.  The  pertinent  event  information  is  given  in  Table  II- 1.  The  first  event  is  one  Of  the 
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TABLE  11-1 

PDE  INFORMATION  FOR  THE  TWO  EVENTS  USED  IN  THIS  STUDY 

Station 

Date 

(1976) 

Origin  Time 

Latitude 

(°N) 

Longitude 

Depth 

Region 

mb 

ANMO 

MAIO 

12  May 
16  March 

19:50:00.2 
07:28:57. 6 

37.2 

27.3 

1 16.  2°W 
55. 1°E 

0 

33 

S.  Nevada 
S.  Iran 

4.9 

5.4 

few  NTS  explosions  which  did  not  overdrive  the  Albuquerque  SRO.  Rotated  LP  seismograms 
are  shown  in  Figs.  II -6  and  II -7,  respectively.  A close  examination  of  Fig.  II -7  reveals  the 
fundamental  and  first  higher  Rayleigh-wave  modes  interfering  on  the  LPZ  component  and  the 
barest  hint  of  a higher  Love-wave  mode  preceding  the  fundamental  on  the  LPT  component. 

In  processing  these  and  other  similar  data,  we  found  that  the  filter  learning  constant  sug- 

7 

gested  by  Griffiths  (a  = 0.75)  was  optimum;  it  gave  the  best  balance  between  slow  adaptation 
on  the  one  hand  and  noisy  performance  on  the  other.  Best  results  were  obtained  by  using  a 
suite  of  filter  lengths  between  8 and  20  sec  long  where  the  data  were  sampled  at  1 Hz.  As  a 
general  rule,  short  filter  lengths  produced  narrower  contours  but  were  less  sensitive  to  longer 
periods  than  the  long  filters. 

The  results  are  presented  in  Figs.  II -8  through  11-10,  which  are  group-velocity-dispersion 

curves  plotted  to  overlie  the  spectral -contour  plots.  Doubtful  values  are  indicated  by  dashed 

lines.  No  attempt  was  made  to  remove  either  the  SRO  instrument  or  source-group  delays. 

Since  maximum  values  of  this  combined  effect  are  less  than  10  sec  for  the  band  indicated  in  the 

plots,  the  group-delay  error  cannot  be  more  than  0.01  km/sec.  Both  events  were  approximately 

1000  km  away  from  the  respective  SROs.  Another  possible  analysis  delay  is  due  to  the  filter 

learning  response  tirpe.  However,  the  examples  run  by  Griffiths  using  the  same  learning  con- 
7 

stant  show  that  this  is  a negligible  effect. 

The  results  in  Figs.  II-8  through  II -10  suggest  that  the  two  propagation  paths,  the  Colorado 
Plateau  in  the  first  case  and  Iranian  Platform  in  the  second,  exhibit  relatively  low  group  veloci- 
ties. Figure  II -8  includes  a comparison  with  the  group -velocity  dispersion  from  Model  72 T 
(see  Ref.  9)  appropriate  to  the  Basin  and  Range  geologic  province.  This  comparison  indicates 

that  the  Colorado  Plateau  structure  is  significantly  different  than  that  for  Model  72 T.  Other 
10 

investigators,  using  WWSSN  LP  recordings  of  the  Rio  Blanco  explosion,  have  found  that  the 

Colorado  Plateau  has  a relatively  thick  crust  (~40  to  45  km)  for  the  Western  U.S.  underlain  by 

the  usual  low -velocity  upper  mantle.  The  group -velocity-dispersion  curve  in  Fig.  II -8  is  quite 

10 

similar  to  one  presented  by  Keller  et  al.  for  the  northern  Colorado  Plateau.  It  therefore  lends 
additional  support  for  their  conclusions. 

Figures  II -9  and  11-10  are  Rayleigh-  and  Love-wave  group-velocity-dispersion  curves  for 
the  Iranian  Platform.  In  this  case,  it  was  possible  to  measure  a Rayleigh-wave  higher  mode 
and  see,  at  least,  the  indication  of  a Love -wave  higher  mode.  The  high-resolution  method  was 
used  to  advantage  because  the  higher  modes  were  intermingled  with  the  fundamental  modes. 

Both  fundamental  modes  exhibit  extraordinarily  low  group  velocities  which  can  be  compared 
with  the  dispersion  from  the  Gutenberg  continent  model  also  included  in  the  figures.  Such  group 
velocities  would  indicate  an  anomalously  thick  crust,  possibly  50  to  60  km.  This  interpretation 


24 


is  consistent  with  the  tectonic  explanation  that  Iran  is  being  compressed  by  the  northward  move- 
ment  of  the  Arabian  lithospheric  plate.  D w>  McCowan 


D.  A RAYLEIGH-WAVE  STRUCTURE  FOR  NOVAYA  ZEMLYA 


11 


The  sequence  of  large  (presumed)  explosions  detonated  in  Novaya  Zemlya**  general 

usually  good  Rayleigh-wave  seismograms.  A pair  of  these  events,  one  from  the  northe 

the  other  from  the  southern  test  areas  on  Novaya  Zemlya,  recorded  at  ALPA  were  anal^ 

12 

using  the  single-station,  two-event  method  due  to  Alexander.  This  ’analysis  produced 


velocity -dispersion  curve  appropriate  to  the  source  region.  The  pertinent  information 
two  events  is  given  in  Table  II -2. 
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TABLE  11-2 

PDE  INFORMATION  FOR  THE  TWO  NOVAYA  ZEMLYA  EVENTS 
USED  IN  THIS  STUDY 

Event 

Date 

Origin  Time 

mb 

Latitudet 

(°N) 

Longitude 

(°E) 

Northern 

28  August  1972 

05:59:57 

6.3 

73.3 

55.1 

Southern 

2 November  1974 

04:59:56 

6.7 

70.8 

54.1 

t These  epicenters  imply  an  event  separation  of  285  km. 

j 

The  method  consists  of  crosscorrelating  the  set  of  seismograms  from  one  event  with  those  ' 
from  the  other.  Since  these  were  array  observations,  the  resulting  set  of  crosscorrelations 
from  corresponding  array  elements  was  then  averaged.  If  the  following  assumptions  arq  valid: 

(1)  The  source  region  can  be  characterized  by  an  average  plane  layered 
structure; 

(2)  The  Rayleigh  waves  from  both  events  traverse  a common  path  from  the 
nearer  (northern)  event  to  the  receiver; 

(3)  Interference  due  to  multipathing  can  be  neglected; 

(4)  The  source  phase  for  the  two  events  is  either  identical  or,  at  least, 
known;  . 


then  the  phase  of  the  resulting  cross  spectrum  will  exhibit  the  phase-velocity  dispersion 
source  region,  i.e.. 


in  the 


exp  [— i 


C(<o)t  • 


Here,  D is  the  separation  between  the  two  events. 

Assumption  (2)  above  is  fortuitously  met  by  this  combination  of  events  and  station.  The 
two  events  are  almost  perfectly  lined  up  on  a great-circle  path  which  traverses  the  Arctic  Ocean 
and  passes  through  ALPA.  Because  this  path  consists  largely  of  an  ocean  basin,  multipathing 
effects  can  be  expected  to  be  small.  However,  as  will  be  seen  in  the  results,  there  were  minima 


TABLE  11-3 

LAYER  PARAMETERS  FOR  MODEL  NZ1 

Thickness 

P-Wave  Velocity 

S-Wave  Velocity 

Density 

2.0 

4.  95 

2.26 

2.65 

3.0 

5.55 

2.73 

2.70 

5.0 

6.  00 

3.03 

2.75 

5.0 

6.00 

3.05 

2.75 

5.0 

6.20 

3.35 

2.  80 

10.0 

6.78 

3.97 

2.85 

10.0 

6.78 

3.97 

2.85 

10.0 

6.78 

3. <89 

2.85 

10.0 

7. 99 

4.33 

3.25 

15.0 

8.10 

4.29 

3.25 

20.0 

8.  10 

4.27 

3.25 

20.0 

8.20 

4.39 

3.28 

40.0 

8.20 

4.52 

3.30 

60.0 

8.20 

4.59 

3.31 

100.0 

8.20 

4.60 

3.32 

0 

8. 40 

4.67 

3.37 

TABLE  11-4 

DATA  AND  COMPUTED  VALUES  FROM  MODEL  NZ1 

Period 

Observed  C(T) 

Computed  C(T) 

Computed  U(T) 

51.08 

3.83 

3.85 

3.45 

42.70 

3.75 

3.77 

3.35 

36.49 

3.69 

3.69 

3.24 

31.93 

3.63 

3.62 

3.15 

28.38 

3.56 

3.56 

3.07 

25.54 

3.50 

3.49 

2.99 

23.22 

3.45 

3.44 

2.92 

21.28 

3.39 

3.38 

2.84 

18. 24 

3.28 

3.28 

2.70 

15.96 

3.  18 

3. 18 

2.57 

14.19 

3.09 

3.09 

2.48 
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in  the  cross -power  spectrum  which  may  have  been  caused  by  multipathing.  Since  there  is  prac- 
tically no  available  information  in  the  open  literature  on  the  crustal  structure  of  Novaya  Zemlya, 
it  is  impossible  to  know  a priori  whether  or  not  assumption  (1)  is  met.  Instead,  the  spirit  of 
this  experiment  is  to  postulate  it  as  being  true,  and  see  whether  or  not  the  results  are  reason- 
able. We  have  also  assumed  that,  both  events  being  presumed  explosions,  they  will  have  iden- 
tical source  phases.  However,  the  possibility  exists  that  either  event  could  have  triggered 
tectonic  stress  release  which  would  add  an  earthquake -like  source  phase  to  the  spectra. 

The  results  are  shown  in  Fig.  II- 1 1 . The  data  points  are  plotted  as  small  dots  with  a 
smooth  curve  drawn  through  them  in  the  lower  graph.  The  upper  graph  is  a plot  of  the  relative 
power  in  the  averaged  cross  spectrum.  As  can  be  seen,  there  is  a pronounced  minimum  near 
the  22 -sec  period;  this  corresponds  to  a dip  in  the  measured  phase  velocities.  Two  other  areas 
of  erratic  phase  velocity  are  at  the  17-  and  30-sec  periods,  respectively.  The  former  of  these 
corresponds  to  another  minimum  in  the  relative  power  of  the  latter  to  a hump-like  sidt  lobe 
above  the  main  peak.  We  did  not  consider  any  of  these  phase -velocity  fluctuations  to  be 
significant. 

The  X's  in  the  lower  graph  are  computed  phase  velocities  from  the  model  given  iq  Ta- 
ble II-3.  The  data,  as  well  as  the  computed  phase  and  group  velocities,  are  given  in  Table  II -4 
where  one  can  see  that  the  rms  error  of  the  fit  to  the  phase -velocity  data  is  approximately 

0.01  km/sec.  This  model  was  obtained  by  using  the  stochastic  inverse  method  applied  to 

13 

surface -wave -dispersion  data. 

The  model  is  characterized  by  a 50-km-thick  crust  and  a slight  shear-wave  low-velocity 
zone  in  the  lower  mantle.  The  latter  is  probably  an  artifact  of  the  inversion  process,  since 
lower  bounds  on  the  standard  deviations  of  the  model  parameters  exceed  this  dip.  The  50-km- 
thick  crust  is  the  result  of  several  attempts  at  data  inversion.  It  minimizes  extreme  oscilla- 
tions in  adjacent  layer  parameters.  This  model  is  not  offered  as  a unique  solution,  as  it  is 
well  known  that  surface -dispersion  data,  band-limited  and  for  only  the  fundamental  mode,  can- 
not determine  a unique  structure.  Furthermore,  differential  model  parameter  changes  within 
the  averaging  kernels  produced  by  the  stochastic  inversion  method  agree  equally  well  with  the 
data  within  its  observational  errors.  We  only  offer  this  as  a plausible  model  that  fits  our  lim- 
ited data. 

The  50-km-thick  crust  of  Model  NZ1  is  consistent  with  the  explanation  that  Novaya  Zemlya 

14 

is  a northern  extension  of  the  Ural  mountains.  Kosminskaya  et  al.  present  results  for  Russian 
crustal  structures  determined  by  "Deep  Seismic  Sounding"  (DSS).  In  particular,  their  lpg.  11 
is  a contour  map  of  Moho  depth  for  the  whole  Soviet  Union.  The  Ural  mountain  belt,  which  gen- 
erally lies  along  the  60° E longitude  meridian,  is  pictured  as  having  a 45-km -thick  crust.  Un- 
fortunately, none  of  their  profiles  are  appropriate  to  Novaya  Zemlya;  however,  the  P-wave 
velocity  structure  of  our  model  is  similar  to  that  for  their  Ukranian  Shield  profile. 

D.  W.  McCowan 
P.  Glovert 
S.  S.  Alexander^ 


t Pennsylvania  State  University,  University  Park,  PA. 
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E.  PURE  PATH  DISPERSION  OF  OVERTONES  OF  SURFACE  WAVES 


Most  of  the  inferences  with  respect  to  the  differences  between  the  continental  and  oceanic 

15-17 

mantle  have  been  made  through  interpretation  of  "pure  path'1  dispersion  of  mantle  waves. 

As  the  resolving  power  of  these  data  becomes  poor  for  depths  greater  than  300  km,  the  details 
of  the  models  below  this  level  are  rather  arbitrary  and  it  has  been  shown  that  they  are  unnec- 
essary to  explain  the  observations.  One  of  the  ways  to  increase  the  overall  resolution  of  the 
upper-mantle  structure  would  be  to  obtain  reliable  measurements  of  dispersion  of  overtones 
to  periods  as  short  as  20  sec. 

Figure  11-12  is  meant  to  illustrate  this  point.  Models  1066A  and  1066B  of  Gilbert  and 
1 8 

Dziewonski  give  an  equally  good  fit  to  the  normal-mode  data,  even  though  their  upper-mantle 

structures  differ  significantly  in  detail;  model  B has  abrupt  discontinuities,  while  model  A is 

smooth.  Group  velocities  for  the  fundamental  mode  are  practically  identical  for  periods 

greater  than  100  sec;  for  shorter  periods,  the  discrepancies  are  caused  by  differences  in  the 

details  of  crustal  structures.  Yet  the  dispersion  curves  for  the  fourth  and  sixth  overtones 

are  significantly  (5  to  7 percent)  different  in  the  range  of  phase  velocities  that  correspond  to  the 

turning  point  between  300  and  700  km  depth.  It  is  clear  that  if  precise  group  velocity  data  were 

available,  much  more  could  be  said  about  the  actual  structure  of  the  upper  mantle.  At  the  same 

time.  Fig.  11-12  illustrates  the  fact  that  such  measurements  could  not  be  made  using  individual 

recordings,  as  there  are  multiple  intersections  of  group  velocity  curves  where  dispersion 

could  not  be  resolved  by  the  frequency-time  analysis. 

Initial  progress  in  the  development  of  techniques  of  mode  separation  has  been  made  by 
19  20  21 

Forsyth  and  Nolet  ' who  have  used  the  " stacking*'  approach  in  order  to  achieve  mode 
separation.  The  procedure  outlined  here  is  different  in  several  aspects  from  those  applied 
previously. 

(1)  Having  a regional  network  of  seismographs  (for  example,  WWSSN  sta- 
tions in  Western  Europe)  and  a given  location  of  the  source  and  its  source 
mechanism,  we  compute  synthetic  seismograms  for  the  first  orbit  (Rl) 
for  each  overtone  whose  dispersion  is  to  be  measured;  this  should  be 

done  for  each  of  the  stations  in  the  array;  the  main  elements  of  the  theory 

2 2 th 

have  been  described  by  Gilbert.  The  synthetic  spectrum  of  the  n over- 
tone is 
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where  r is  the  position  of  the  receiver  in  the  epicentral  coordinate  sys- 
tem, rn  is  the  radius  of  the  source,  a is  the  Earth's  radius,  a(a>)  is  the 
0 n 

attenuation  factor  (a  = a>/2Q),  A is  the  epicentral  distance,  nu(u>)  is 

the  group  velocity,  A.  (r,  rn,co)  are  analogous  to  those  in  Gilbert  and 
18  n l ^ u 

Dziewonski's  expressions  (2.1.30)  for  spheroidal  modes  and  (2.1.31) 

for  toroidal  modes,  with  the  functions  xfn(9)  replaced  by  those  in  Eq.  (21) 
22  1 

of  Gilbert;  T(oo)  are  the  spectra  of  the  six  independent  components  of 
the  moment  rate  tensor. 
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(2)  Compute  the  stack  of  crosscorrelograms  between  the  actual  recordings 
(Sj.  where  j is  the  index  of  a station  within  the  array)  and  the  synthetic 
seismogram  for  a given  overtone  (nwj)-  The  spectrum  of  such  a stack  is 

nr(«)  = Tj  Sj(w)  nw.(w) 
j 

where  the  bar  denotes  the  complex  conjugate.  This  process  can  be  re- 
peated for  any  number  of  events  from  roughly  the  same  source  region. 

(3)  Measure  differential  group-  and  phase -velocity  dispersion  using  the 
"residual  dispersion  method"  of  Dziewonski  et  al.^ 

Figures  II-13(a-p)  show  results  of  a synthetic  experiment  in  which  an  attempt  was  ifnade  to 
isolate  the  energy  associated  with  individual  overtones  using  an  earthquake  in  Japan  (49l05°N, 
153. 6°E,  depth  135  km)  and  15  WWSSN  stations:  AQU.  ATU,  COP.  ESK,  HLW,  1ST,  JER.  KEV, 
KON,  MAL,  NUR,  PTO,  STU,  TOL,  and  UPP.  The  observed  seismograms  were  slmi|lated 
by  superposition  of  fundamental  Rayleigh  and  Love  modes  and  their  first  seven  overtones. 

The  period  range  of  the  analysis  extended  from  20  to  500  sec.  The  retrieved  signal  for  the 
fundamental  Rayleigh  mode  is  nearly  perfect;  in  all  remaining  cases,  we  note  various  degrees 
of  interference  with  other  modes.  By  proper  windowing  and  truncation,  it  should  be  possible 
to  recover  most  of  the  information  for  seems  to  be  mostly  disturbed  by  a high- 

frequency  energy  (the  third  overtone  is  continuous  with  the  second  after  intersection  wi  ;h  the 
core-mantle  boundary  Stoneley  branch).  The  same  is  true  of  /^S^;  is  quite  good, 

but  and  might  not  yield  reliable  results.  Of  the  toroidal  modes,  probably 

^Tj  has  the  best  signal-to-noise  ratio,  but  several  others  should  also  prove  useful. 

We  would  like  to  point  out  that  these  results  are  for  a single  source.  If  several  sources 
with  slightly  different  locations  were  used,  the  phase -equalization  procedure  would  always 
lead  to  reinforcement  of  the  signal  of  interest,  while  other  modes  having  different  phasfe-arrival 
times  would  tend  to  cancel. 

The  method  can  be  applied  to  investigate  purely  continental  paths  under  the  Eurasia  and, 
for  example,  oceanic  paths  under  the  Pacific  using  an  array  of  stations  in  western  United  States 
and  Canada.  Also,  stations  in  the  eastern  United  States  and  Canada  could  be  used  to  study  dis- 
persion from  suitably  located  events  on  the  Mid-Atlantic  Ridge,  and  European  stations  Would 
have  a purely  oceanic  path  for  the  events  in  the  Caribbean. 

If  these  studies  are  successful,  and  there  is  every  indication  that  they  should  be,  it  would 
be  very  desirable  to  apply  a relocatable  network  of  broadband  digital  instruments  for  tl}e  spe- 
cific purpose  of  investigating  the  multimode  surface -wave  dispersion  along  selected  paljhs. 

T.  A.  Chout 
A.  M.  Dziewonski 


t Department  of  Geological  Sciences,  Harvard  University,  Cambridge,  MA  02138. 
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Fig.  II-l.  Frequency  shift  by  -u>c  to  achieve 
bandwidth  reduction. 
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Fig.  II -2.  Amplitude  and  phase  characteristics  of  filtered  and  frequency-shifted 
waveform  resulting  from  analysis  of  oceanic  surface-wave  data  for  center  fre- 
quency = 0.04  Hz  (T  = 25  sec). 
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Fig.  II-3.  Group  and  phase  velocities  for  well-dispersed  continental  seismic 
data  determined  using  modified  Dziewonski  method.  (Correction  for  source 
phase  not  included. ) 
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Fig.  II -4.  Polar  plot  centered  at  Mashad  showing  events  recorded  at 
Mashad  SRO  to  which  group  velocity  has  been  determined.  Numbers 
refer  to  events  from  which  group  velocity  is  shown  in  Fig.  II -5. 
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Fig.  II — 5.  Group -’velocity  dispersion  measured  at  Mashad  for  paths 
from  events  numbered  in  Fig.  II-4. 
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Fig.  II-6.  Rotated  LP  seismograms  for  12  May  1976  NTS  explosion 
recorded  at  Albuquerque  SRO. 
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Fig.  II-7.  Rotated  LP  seismograms  for  16  March  1976  Southern  Iran 
earthquake  recorded  at  Mashad  SRO. 
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Fig.  II-8.  Rayleigh-wave  group -velocity 
dispersion  for  12  May  1976  NTS  explo- 
sion. Also  included  for  comparison  is 
dispersion  for  Model  72 T appropriate  to 
Basin  and  Range  structure. 


Fig.  II-9.  Rayleigh-wave  group-velocity 
dispersion  for  16  March  1976  Southern 
Iran  earthquake.  Also  included  for  com- 
parison is  dispersion  for  Gutenberg  con- 
tinent model. 
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Fig.  n-10.  Love -wave  group- velocity 
dispersion  for  16  March  1976  Southern 
Iran  earthquake. 


Fig.  11-11.  Phase-velocity-dispersk' 
suits  for  Novaya  Zemlya.  Lower 
shows  observed  data  (dots),  smooth^ 
terpretation  (line),  and  results  corr 
from  Model  NZ1  (X’s).  Upper  gra; 
relative  power  in  cross  spectrum. 
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Fig.  H-12.  Comparison  of  group-velocity- 
dispersion  curves  for  several  overtones  of 
Love  waves  computed  for  earth  models 
1066A  and  1066B  (Ref.  18). 
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Fig.  II-13(a-p).  Synthetic  examples  of  an  attempt  at  extraction  of  residual  dispersion 
for  16  overtones  of  surface  waves  by  processing  synthetic  3 -component  records  from 
15  European  WWSSN  stations. 
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III.  MAGNITUDE,  YIELD,  AND  ENERGY 


A.  AZIMUTHAL  P-WAVE  AMPLITUDE  VARIATIONS  FROM  LARGE  EXPLOSIONS 

In  recent  years,  increasing  use  has  been  made  of  body  phase  amplitudes  to  constra|n  source 
and  structure  models.  These  constraints  often  are  applied  with  the  barest  understanding  of  am- 
plitude variations  from  the  effects  of  anelasticity  and  scattering.  These  effects  are  difficult  to 
measure  independently;  at  best,  both  effects  are  lumped  together  by  assigning  an  effective  Q 
or  t*  to  the  particular  ray  path. 

Teleseismic  signals  in  the  band  pass  of  SP  instruments  vary  by  a factor  of  4 globally  due  to 
the  effects  of  station  bias  (North  ).  In  some  regions,  these  biases  change  dramatically  over  dis- 
tances that  are  short  on  a teleseismic  scale,  e.g.,  the  station  BKS  at  Berkeley,  California  and 
TFO  in  Arizona  which  have  an  arc  separation  of  about  6°  and  have  biases  of  +0.18  and  —0.33, 
respectively,  in  units  of  log  (A/T).  The  equivalent  amplitude  differences  are  about  a factor  of 
3 for  signals  of  about  1-sec  period.  Presumably  similar  differences  would  occur  over  shorter 
distances  if  stations  were  more  densely  distributed. 


1 


North  showed  that  biases  and  the  regional  tectonic  setting  are  correlated,  but 
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there  is  little  azimuthal  variation  in  these  biases  for  a given  region  or  for  a given  station.  Yet 

azimuthal  variation  related  to  near-source  effects  is  well-known  from  studies  of  SP  signals 

2 3 

from  explosions  and  from  earthquakes,  particularly  those  in  island  arc  regions. 

We  have  begun  a study  of  the  amplitude  variation  of  P waves  as  a function  of  azimuth  with 
LP  amplitudes  from  two  explosions:  Cannikin  in  Amchitka  Island  [Fig.  III-l(a)],  and  on£  in 
Novaya  Zemlya  [Fig.  III-l(b)].  The  LP  amplitudes  are  particularly  useful  for  this  study  because 
the  signal  is  simple,  e.g.,  see  Frasier,  Sec.  H below.  As  a result,  the  chance  of  reading  errors 
is  minimized  and  the  formal  modeling  of  these  waveforms  is  encouraged.  The  object  of  this 
work  is  to  separate  the  gross  effects  of  scattering,  including  focusing  and  defocusing,  fjrom 
those  of  anelasticity  and  thus  to  measure  the  real  t*  for  the  earth.  Explosions  are  preferred 
to  earthquakes  for  this  work  because,  ideally,  explosions  radiate  seismic  energy  uniformly, 
although  some  deviation  from  uniformity  may  result  from  the  release  of  tectonic  strain 

Figure  III- 1 (a)  shows  that  P waves  from  Cannikin  arrive  at  European  and  extreme  eastern 
North  American  stations  with  amplitudes  diminished  by  a factor  of  2 to  3 from  the  average 
amplitude  recorded  elsewhere  (about  7 normalized  to  an  arc  distance  of  60°).  This  afimuthal 
variation  can  be  explained  by  defocusing  of  signals  that  pass  through  the  seismic  zone  ihclined 

4 

from  the  Aleutian  Arc.  Three-dimensional  ray  tracing  by  Jacob  has  shown  that  body  phases 
recorded  at  European  stations  from  the  Amchitka  explosion,  Longshot,  will  travel  through  the 
slab  and  show  negative  travel-time  residuals.  We  intend  to  quantify  this  amplitude  effect  in 
the  future. 

The  smaller-magnitude  Novaya  Zemlya  explosion  shows  no  broad  amplitude  variation  com- 
parable with  the  Cannikin  signals;  however,  amplitudes  recorded  in  the  Middle  East  ovier  an 
azimuth  window  of  about  20°  are  larger  than  a mean  amplitude  for  this  event  by  about  a factor 
of  2.  In  the  middle  of  this  window  are  two  stations,  AAE  and  NAI,  near  the  East  African  rift 
that  recorded  low  amplitudes.  The  azimuth  window  corresponding  to  stations  in  the  Western 
United  States  as  far  east  as  Texas  is  also  characterized  by  low  amplitudes.  These  observations 
are  consistent  with  North’s  study  in  that  stations  near  the  rift  and  those  in  the  Western  United 
States  show  m^  biases  equivalent  to  amplitude  diminutions  of  60  and  40  percent,  respectively. 

T.  J.  Fitch 
C.  W.  Frasier 


B.  SEISMIC  SCALING  INFORMATION  FROM  SRO  DATA 


Seismic  Scaling  Laws  are  of  great  importance  to  discrimination  since  determination  of  the 
correct  law  would  allow  one  to  predict  M -m,  , M -Moment,  and  even  m,  -Moment  relations  for 

SOS  £ D 

earthquakes  and  explosions.  Unfortunately,  as  has  been  pointed  out,  measurements  of  magni- 
tude and  moment  are  often  insufficiently  precise,  primarily  because  of  our  lack  of  knowledge 
about  source  mechanism  (as  opposed  to  size)  and  propagation  effects,  to  enable  the  various  laws 
proposed  to  be  tested.  A powerful  way  of  testing  these  laws  is  to  compare  spectral  ratios  of 
events  of  different  sizes,  but  from  the  same  source  region, with  model  predictions.  The  SRO 
stations,  with  their  large  dynamic  range  and  high  digital -data  quality,  are  particularly  well- 
suited  to  the  calculation  of  spectral  ratios,  and  we  describe  here  the  collection  and  preliminary 
analysis  of  a unique  data  set  for  this  purpose. 

On  21  January  1976,  a large  earthquake  occurred  in  the  southern  Kurile  Islands.  During 
the  following  20  days,  some  120  aftershocks  of  this  event  were  located  by  the  USGS.  Both  sur- 
face and  body  waves  from  most  of  these  events  were  well-recorded  at  Mashad  SRO:  the  two 
other  SROs  installed  at  the  time  were  either  out  of  operation  (Albuquerque)  or  of  poor  data  qual- 
ity (Guam).  All  the  aftershocks  reported  have  been  relocated  using  a master-event  technique, 
with  the  main  shock  as  the  master:  use  of  this  method  considerably  reduces  the  aftershock  area 
(Fig.  III-2).  A remarkable  feature  of  this  sequence  is  the  number  of  depth  phases  reported:  for 
over  60  percent  of  the  events,  three  or  more  (pP-P)  times  constrain  the  depth  to  45  to  55  km. 

We  thus  have  a set  of  events  covering  a large  range  of  magnitude,  located  within  a very 
small  area,  and  well-recorded  at  Mashad  SRO.  The  chief  remaining  uncertainty  is  that  of  the 
source  mechanism:  most  of  the  events  were  too  small  for  first  motion  to  be  observable  on  WWSSN 
seismograms.  Those  which  are  identifiable  show  only  compression,  probably  indicating  faulting 
of  thrust  type,  which  is  consistent  with  the  tectonics  of  this  area.  Perhaps  the  most  convincing 
evidence  of  a uniform  fault  plane  orientation  for  these  events  is  the  remarkably  similar  nature 
of  all  the  short-period  P waveforms  recorded  at  Mashad. 

Figure  III— 3 shows  vertical -component  Rayleigh-wave  trains  from  five  of  the  Kurile  events, 
as  recorded  at  Mashad  SRO.  Records  (a)  and  (e)  are  contaminated  by  surface  waves  from  smaller 
events  in  the  same  source  region.  This  figure  dramatically  illustrates  the  dynamic  range  of  the 
SRO  recording  system:  peak-to-peak  Airy  phase  amplitudes  corrected  for  instrument  response, 
range  from  720  p (a)  to  0.36  \i  (e),  corresponding  to  Mg  = 7.6  and  4.2,  respectively.  We  plan 
eventually  to  utilize  sophisticated  filtering  techniques,  such  as  polarization  filtering  (see 
Sec.  IV-D  of  this  report)  and  time  variable  filtering  based  on  the  easily  measurable  group- 
velocity-dispersion  characteristics  of  the  path  to  improve  signal-to-noise  ratios  and  eliminate 
contaminating  wave  trains  from  sources  in  the  same  region  or  elsewhere. 

Figure  III-4  shows  amplitude  spectra  of  events  (a)  through  (d)  of  Fig.  III-3.  The  spectra 
are  calculated  from  exactly  the  same  time  window  in  each  case;  thus,  the  only  difference  be- 
tween them  should  be  caused  by  seismic  scaling  effects.  The  effects  of  the  source  mechanism, 

propagation  path,  and  instrument  should  be  identical.  Events  (a)  through  (d)  are  of  M = 7.6, 

2 6 s 
6.6,  5.6,  and  4.6,  respectively.  The  uj  law  of  Aki  predicts  that  these  Mg  values  correspond 

to  Mq  = 4 x 1027,  3.5  x 102^,  2 X to2"*,  and  1.8  x 102^  dyn-cm  and  corner  periods  t = 105,  30, 

10,  and  2 sec,  respectively.  Thus,  within  the  period  range  over  which  we  have  calculated  the 

spectra,  we  would  expect  to  see  significant  changes  in  spectral  shape  over  the  magnitude  range 

considered. 
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Figure  III-5  shows  spectral  ratios  of  events  (a/b),  (b/c),  and  (a/c)  calculated  at  periods 

of  10,  20,  50  and  100  sec  (crosses).  The  solid  line  shown  is  the  spectral  ratio  predicted  by  the 
2 

cj  model.  At  the  shorter  periods  of  10  and  20  sec,  the  agreement  between  observation  and 
theory  is  quite  good,  but  both  the  ratios  involving  event  (a)  are  considerably  less  than  predicted 
at  the  longer  periods. 


These  data  do  not  therefore  satisfy  the  <x>  law:  we  intend  to  investigate  other  scaling  laws 
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h as  those  involving  ar 
quality  of  the  spectra  obtained. 
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(such  as  those  involving  and  falloff  above  the  corner  frequency)  and  to  improve  the 


R.  G.  North 
C.  W.  F rasier 


C.  A RAPID  AND  INEXPENSIVE  PROCEDURE  FOR  MEASURING 

THE  RESPONSE  OF  THE  WWSSN  SP  INSTRUMENT 

As  installed,  the  WWSSN  SP  recording  system  provides  for  a "weight  lift"  calibration  pulse. 
Unfortunately,  at  the  operating  drum  speed,  the  pulse  cannot  be  resolved  in  time.  In  fjact,  it  is 
difficult  to  see  at  all.  Consequently,  the  preferred  calibration  procedure  recommended  by  the 
manufacturer  (Geotech)  is  to  drive  the  calibration  coil  with  AC  currents  and  measure  the  am- 
plitude and  phase  delay  of  the  resulting  output  relative  to  the  input  current.  To  do  this  properly 
requires  electronic  equipment.  Alternatively,  the  seismometer  can  be  removed  and  nr  ounted 
on  a shake  table  where  its  amplitude  and  phase  responses  can  be  measured  relative  to  Ithe  de- 
flection of  the  table.  While  this  latter  procedure  eliminates  all  seismic  noise  from  the  response 
measurement,  it  too  requires  specialized  equipment.  By  far,  the  simplest  procedure  {is  to 
speed  up  the  drum  to  resolve  the  calibration  pulse.  Then  it  can  be  Fourier  analyzed  in  much 
the  same  manner  as  are  the  LP  calibration  pulses.  We  present  here  a relatively  simpjle  means 
of  doing  this  with  a geared-down  electric  motor. 

The  motor  and  pulley  arrangement  is  shown  in  Fig.  HI-6.  The  Bodine  motor  has  a.  10-rpm 
reduction  drive  which  generates  72  in.-oz  of  torque.  In  our  experiment,  the  high  torque  was  a 
useful  feature  because  of  friction  in  the  recording  drum  drive  system.  To  this  motor  ^vas  added 
a 2 -in. -diameter  pulley  with  an  O ring  to  improve  traction.  The  resulting  recording  ifate  was 
therefore  approximately  1 in. /sec. 

We  are  indebted  to  the  Weston  Seismic  Observatory  (particularly  Mr.  R.  O.  Ahner)  for  allow- 
ing us  to  try  out  our  procedure  on  their  SPZ  recorder  during  a record -changing  sessicjn.  In 
all,  we  recorded  11  calibration  pulses  with  roughly  an  even  mix  of  both  polarities,  of  ^hich  10 
were  usable.  These  were  then  hand -digitized  on  our  Bendix  digitizer  at  approximately  28-Hz 
sampling  rate.  The  10  sets  of  digital  data  were  then  phaseless  low-pass  filtered  and  resampled 
at  20  Hz.  Each  trace  was  aligned  by  eye  and  summed  to  produce  the  pulse  shown  in  Fig.  III-7. 
The  measured  overshoot  ratio  is  approximately  12  to  1. 

The  seismometer  used  in  the  WWSSN  SP  system  is  a variable  reluctance  Benioff  ihstrument. 

7 

As  such,  the  well-known  Hagiwara  formulation  of  the  system  response  is  inappropriate  because 
the  inductance  of  the  seismometer  has  been  neglected.  The  correct  form  of  the  equatipns  has 
been  derived  by  Chakrabarty  and  Choudhury8  and,  when  applied  to  the  WWSSN  SP  system,  pre- 
dicts noticeably  different  amplitude  and  phase  responses  than  zero-inductance  approximations. 

In  particular,  the  phase  varies  over  5tt/2  rad  from  DC  to  infinite  frequency  as  opposed  to  2ir  rad 
in  zero-inductance  systems.^  The  installed  overshoot  ratio  for  the  WWSSN  SP  system  is 
supposed  to  be  17  to  1. 
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Amplitude  and  phase  spectra  computed  from  the  calibration  pulse  in  Fig.  Ill- 7 are  shown 
in  Figs.  Ill— 8 and  III-9,  respectively.  The  amplitude  response  is  usable,  with  some  smoothing, 
from  DC  to  5 Hz  but  the  phase  response  is  well-determined  only  between  DC  and  3 Hz.  The 
discrepancy  in  usable  bandwidths  is  explainable  in  terms  of  the  fiducial  phase  delay  of  the  pulse. 
Since  the  sampling  frequency  is  20  Hz,  the  phase  error  due  to  a starting-time  error  of  l/40  sec 
(half  the  sampling  interval)  at  3 Hz  is  approximately  0.5  rad.  The  amplitude  response  is  inde- 
pendent of  this  effect.  The  DC  phase-angle  asymptote  of  —3 tt/ 2 rad  predicted  by  the  theory  is 
supported  by  the  results  in  Fig.  III-9 • Unfortunately,  due  to  the  limited  bandwidth  of  the  experi- 
ment, it  is  not  possible  to  verify  a high-frequency  asymptote  of  +7r  rad.  Judging  from  the  mea- 
sured overshoot  ratio,  this  particular  system  appears  to  be  underdamped  relative  to  the  speci- 
fication. Whether  or  not  the  actual  system  parameters  can  be  obtained  through  analysis  of  these 

10 

measurements,  as  Mitchell  and  Landisman  have  done  for  WWSSN  LP  systems,  remains  to  be 
seen.  We  have  merely  demonstrated  that  it  is  feasible  to  run  the  recording  drum  faster  and 
record  accurate  calibration  pulses.  D w>  McCowan 

L.  E.  Sargent 


D.  ESTIMATION  OF  STATION  DETECTION  CHARACTERISTICS 

FROM  BULLETIN  DATA 

1 1 

In  the  previous  SATS,  a study  was  described  in  which  a maximum-likelihood  method  was 
used  to  determine  the  station  detection  parameters  for  a set  of  72  stations,  using  data  from  the 
Bulletin  of  the  International  Seismological  Center  (ISC)  for  the  years  1964  through  1973. 

The  present  investigation  is  a continuation  of  the  earlier  one,  and  is  aimed  at  refining  the 
methods  that  can  be  used  to  extract  detection  parameters  from  these  kind  of  data.  Particular 
emphasis  has  been  placed  on  documenting  temporal  changes  in  the  detection  threshold  of  some 
existing  stations,  on  examining  log  A/T  thresholds  as  compared  with  magnitude  thresholds,  on 
obtaining  further  information  about  the  seismicity  parameter  /3,  and  on  the  inclusion  of  satura- 
tion thresholds. 

First,  a pass  was  made  through  the  10  years  of  ISC  data,  and  counts  of  log  A/T  values  were 
extracted  for  each  of  the  years  for  the  72  stations  listed  by  North  (see  p.  2 in  Ref.  11).  Com- 
parison of  yearly  plots  of  frequency  against  log  A/T  revealed  substantial  variations  of  detection 
characteristics  with  time.  In  a few  cases,  this  was  attributable  to  changes  in  station  magnifi- 
cation. However,  in  most  instances  the  variations  had  no  obvious  explanation,  and  were  pre- 
sumed to  be  due  to  variable  operator  performance.  Only  12  of  the  72  stations  were  consistent 
enough  that  all  10  years  of  data  could  be  used  in  the  determination  of  detection  parameters,  and 
about  20  of  the  stations  showed  such  erratic  behavior  that  they  were  not  usable.  In  the  remain- 
ing cases,  portions  of  the  catalog  were  used  during  which  the  stations  were  relatively  consistent. 
The  lengths  of  these  portions  generally  lie  in  the  range  4 to  6 years. 

Frequency  against  log  A/T  plots  are  shown  for  three  good  stations  in  Figs.  Ill -1 0 (a) , (b), 
and  (c).  A typical  feature  of  all  these  plots  is  an  extremely  linear  segment  which  is  interpreted 
as  being  due  to  background  seismicity.  These  straight  segments  were  extended  in  both  direc- 
tions, and  used  to  estimate  the  50-percent  detection  threshold  and  its  associated  spread  y^, 
and  the  50 -percent  saturation  threshold  Gg  and  its  spread  y^.  These  parameters,  together 
with  the  station  biases  estimated  by  North,11  are  listed  in  Table  III- 1 for  a particular  subset 
of  30  stations. 
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TABLE  lll-l 

DETECTION  AND  SATURATION  PARAMETERS  FOR  A SET  OF  30  STATIONS 

(Biases  from  North 

Station 

Bias 

gd 

rD 

Gs 

ys 

ALQ 

-0  20 

1.02 

0.30 

2.7 

0,3 

BMO 

-0.29 

0.28 

0.30 

3.1 

0,3 

BNS 

+0.20 

1.61 

0.30 

3.4 

0,2 

BOZ 

-0.06 

0.98 

0.33 

2.3 

0.2 

CAN 

-0.02 

1 .57 

0.40 

3.2 

0.2 

CLL 

+0.20 

1.30 

0.25 

3.1 

0.3 

COL 

+0.01 

1.00 

0.30 

3.0 

0.3 

CPO 

-0.07 

0.87 

0.25 

3.4 

0.4 

EDM 

+0.37 

1 .63 

0.20 

3.0 

0.2 

EKA 

0.00 

1.12 

0.40 

2.5 

0.3 

EUR 

-0.24 

0.67 

0.40 

2.3 

0.4 

GIL 

-0.04 

0.97 

0.45 

2.9 

0.2 

HYB 

+0.19 

1 .62 

0.35 

3.3 

0.3 

KHC 

+0.10 

1 .20 

0.25 

2.9 

0.2 

KON 

+0.07 

1.36 

0.30 

2.7 

0.2 

LON 

-0.30 

1.13 

0.30 

2.9 

0.2 

LPS 

+0.04 

1 .20 

0.25 

2.9 

0.4 

MUN 

+0.15 

1.61 

0.30 

3.1 

0.2 

NDI 

+0.33 

1.77 

0.40 

3.0 

0.2 

NEW 

+0.05 

1.23 

0.35 

2.6 

0.2 

NOR 

-0.14 

1 .02 

0.28 

2.9 

0.4 

PMR 

-0.08 

1 .03 

0.33 

3.4 

0.2 

PNT 

+0.13 

1.40 

0.30 

2.8 

0.2 

POO 

+0.17 

1.54 

0.25 

2.8 

0.2 

PRU 

+0.04 

1.25 

0.20 

3.1 

0.2 

SJG 

+0.24 

1.46 

0.35 

3.1 

0.2 

TFO 

-0.32 

0.26 

0.50 

3.0 

0.3 

TUC 

-0.14 

1.21 

0.30 

2.5 

0.2 

UBO 

-0.11 

0.45 

0.50 

3.1 

0.3 

WMO 

-0.17 

0.41 

0.23 

2.7 

0.5 

The  slopes  of  the  straight  portions  of  these  plots  showed  some  significant  variations.  Ap- 
proximately half  the  stations  studied  contained  slopes  that  were  extremely  close  to  1.00.  These 
included  all  stations  in  the  U.S.,  Alaska,  and  Canada,  together  with  most  stations  in  India  and 
Australia.  On  the  other  hand,  with  occasional  isolated  exceptions,  stations  in  Africa  had  slopes 
in  the  range  1.10  to  1.20,  while  stations  in  Scandinavia,  Europe,  and  Greenland  (19  stations) 
had  slopes  in  the  range  1.20  to  1.30.  An  example  of  one  of  these  is  shown  in  Fig.  Ill- 1 0(c). 
These  very  interesting  observations  are  being  investigated  further.  The  stations  listed  in 
Table  III-l  all  were  consistent,  with  a seismicity  slope  of  1.00.  • 

One  possible  explanation  of  these  slope  variations  is  that  they  are  due  to  the  inclusion  of 
all  events  detected  by  the  station,  regardless  of  distance  and  depth.  This  might  result  in  un- 
usual emphasis  of  certain  portions  of  the  local  spatial  distribution  of  seismicity.  To  test  this, 
another  pass  was  made  through  the  ISC  catalog,  this  time  retaining  only  events  in  the  distance 
range  30°  to  90°,  and  recording  the  station  value  (i.e.,  including  distance-depth  corrections). 
The  frequency-station  m^  plots  for  the  same  three  stations  in  Figs.  Ill-lO(a-c)  are  shown  in 
Figs.  III-ll(a),  (b)  and  (c).  Although  there  is  a slight  increase  in  the  scatter  of  the  data  points, 
the  general  agreement  of  these  plots  with  the  log  A/T  plots  in  Fig.  Ill- 1 0 is  excellent.  The 
slopes  of  the  linear  portions  agree  within  the  scatter  of  the  data,  and  detection  and  saturation 
thresholds  are  consistent.  There  is  a suggestion  that  the  detection  curves  for  station  m^  do 
not  fall  off  as  fast,  at  low  magnitudes,  as  those  for  log  A/T.  This  is  unexplained  at  present. 

Preliminary  conclusions  from  this  experiment  are  that  Bulletin  data  can  be  used  directly 
in  the  determination  of  station  operating  characteristics  and  station  detection  and  saturation 
parameters.  No  explanation  has  yet  been  found  for  the  strong  variation  in  the  seismicity  pa- 
rameter observed  in  different  areas.  In  view  of  the  shapes  of  many  of  the  frequency-station 
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may  lead  to  unreliable 


magnitude  plots,  the  maximum-likelihood  method  described  earlier 
results. 


M.  A.  Chinnery 
J.  C.  Johnston 


E.  COMPUTER  SIMULATION  OF  NETWORK  PERFORMANCE 

The  30  stations  which  are  listed  along  with  their  detection  and  saturation  parameters  in 

Table  III- 1 can  be  viewed  as  a model  of  a global  network.  The  catalog  simulation  program  de- 
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scribed  in  the  previous  SATS  has  been  adapted  to  include  saturation  parameters,  and  it  has 
been  used  to  generate  several  artificial  earthquake  catalogs  using  these  parameters.  A value 
of  3.8  has  been  added  to  all  thresholds  to  simulate  a typical  distance-depth  correction:  this 
allows  the  catalogs  to  simulate  magnitude  determination. 

Events  input  to  the  program,  based  on  a pseudorandom  number  generator  subroutine,  are 
constrained  to  have  a probability  distribution  of  the  form  log  N = a + bm,  with  the  parameter 
b = 1.00.  The  actual  distribution  of  these  events  for  a catalog  with  15,000  detected  events  is 
shown  in  Fig.  Ill- 1 2 . Detections  at  each  station  are  generated  according  to  the  detection  prob- 
ability curve  at  the  station,  and  at  least  3 station  detections  are  required  before  a network  de- 
tection is  declared.  The  distribution  of  network-detected  events  is  also  shown  in  Fig.  Ill  - 1 2 as 
a function  of  the  true  (input)  magnitudes.  The  50 -percent  detection  threshold  for  the  network 
is  4.35,  and  its  50-percent  saturation  threshold  is  about  7.07. 

One  interesting  aspect  of  Fig.  Ill- 1 2 is  the  sharpness  of  the  falloffs  in  detections  at  both 
low  and  high  magnitudes.  The  effective  y for  both  detection  and  saturation  is  less  than  0.2, 
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those 


i.e.,  smaller  than  any  of  the  y’s  for  the  individual  stations.  This  is  interpreted  as  being)  a re- 
sult of  the  requirement  that  3 station  detections  are  necessary  before  an  event  is  declared  as 
a network  detection.  This  unexpected  result  is  being  investigated  further. 

The  artificial  catalog  may  also  be  used  to  compare  true  event  magnitudes  with  network 
magnitudes  estimated  in  the  usual  way,  by  simply  averaging  the  station  magnitudes  from 
stations  detecting  the  events,  without  regard  to  station  bias.  Figure  III- 1 3 shows  incremental 
plots  of  frequency  against  true  magnitude  and  network  average  magnitude.  Notice  the  w|de  dif- 
ference in  these  plots  at  low  magnitudes. 

The  problem  of  bias  in  network  average  magnitudes  is  illustrated  more  clearly  by  Fig.  Ill  - 1 4 
which  plots  average  values  of  magnitudes  as  a function  of  true  magnitude.  The  upper  magnitude 
limit  behaves  as  expected;  since  stations  receiving  larger -than-average  station  magnitudes 
from  an  event  saturate,  the  network  average  is  biased  smaller  than  the  true  magnitude.  At  low 
magnitudes,  the  network  average  magnitude  is  expected  to  be  larger  than  the  true  magnitude. 
However,  the  opposite  is  observed,  and  this  is  attributed  to,  the  negative  biases  of  the  stations 
with  the  lowest  detection  thresholds  (namely  the  VELA  arrays  BMO,  TFO,  UBO,  and  WMO, 
all  in  the  Western  U.S.). 

The  latter  observation  is  likely  to  apply  to  both  the  ISC  catalog  and  the  USGS  catalog  of 
global  earthquakes,  at  least  for  events  within  the  teleseismic  range  (30°  to  90°)  from  thje  West- 
ern U.S.,  which  includes  much  of  the  Pacific  earthquake  belt.  These  same  stations  have*  the  low- 
est detection  threshold  of  all  stations  reporting  to  these  catalogs.  We  therefore  infer  that  the  ap- 
parent detection  threshold  of  these  global  networks  may  be  several  tenths  of  a magnitude  unit 
too  low.  Further  applications  of  computer  simulations  to  the  study  of  global  network  pejr 
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mance  are  continuing. 


M.  A.  Chinnery 


F.  AN  APPLICATION  OF  THE  ENERGY-MOMENT  TENSOR  RELATION 

TO  ESTIMATION  OF  SEISMIC  ENERGY  RADIATED  BY  LINE  SOURCES 

1 1 

In  a previous  SATS,  we  computed  total  radiated  seismic  energies  for  point- source  models 


of  the  I960  Chilean  and  1964  Alaskan  earthquakes.  The  expressions  derived  for  the  cas 
three  different  source  time  functions  are: 
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Here  we  show  how  to  extend  these  results  to  a source  model  consisting  of  a finite  propagating 
line  source. 

Dziewonski  and  Romanowicz  (Sec.  IV-E  of  this  report)  derived  expressions  for  excitation  of 
normal  modes  by  a line  source  that  may  be  approximated  by  a fragment  of  a great  circle,  of 
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angular  extent  4>,  connecting  the  beginning  and  the  end  points  of  the  fault.  They  employed  a 

specific  rotation  of  the  spherical  coordinate  system,  first  suggested  in  a different  context  by 
1 2 

Backus.  Upon  rotation,  the  fault  represents  a fraction  of  the  equator  where  the  origin  point 
has  a colatitude  n/2,  longitude  0,  and  the  end  point  is  at  7r/24>.  Assuming  that  the  source 
mechanism  does  not  change  during  the  process  of  rupture,  the  expression  for  the  spectral  am- 
plitude of  a particular  normal  mode  k is 


+1 


= E §km(^  ^km(u,) 


*kWk  m=-i 


where  C^(cj)  is  the  resonance  curve  and 

J £V>  = Mo(W):Skm(ro>  ir/2)  ■ 
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If  the  moment  release  along  the  source  line  changes  as  (tt/24)  sin  (7r<p/4>),  then  the  excita- 
tion function  is 


* “(«>  = Mo(u,):Skm(ro,  n/Z)  • 
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7T  COS  X 
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M (co)  is  the  spectrum  of  the  moment  rate  tensor  of  the  equivalent  point  source.  Most  fre- 

~ 0 13 

quently  it  is  represented  by  a constant  Mq(o;)  = Mq  (cf.  Kanamori  ),  but  it  could  be  modeled  to 

include  the  effect  of  the  finite  fault  width,  which  would  lead  to  an  asymptotic  2 decrease  in 

1 4 

amplitudes  at  high  frequencies  (Aki  ). 

S^fr  , 7r/2)  is  the  strain  tensor  at  coordinates  (r  tt/2,  0),  where  r is  the  constant  source 
radius  (depth).  The  expressions  for  = e^ir^  71 /^)  can  be  found  in  Gilbert  and  Dziewonski 
[Eq.(2.1.19)  for  the  spheroidal  modes, and  Eq.(2.1.20)  for  the  toroidal  modes]. 

Note  that,  for  (p  = 0, 


Y£m(e,  <p)  = X(m(Q) 


,rn  r2l_+J  ( t --  m)'  1/2  m 
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computations  of  X^m  for  0 = 7r/2  are  particularly  efficient. 
The  term 
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where 


(III-6) 


is  analogous  to  the  directivity  term  of  Ben-Menahem^  [Eq.(2-19)].  Dziewonski  and  Romanowicz 
discuss  this  analogy  in  some  detail  in  Sec.  IV-E  of  this  report. 

XL 

Thus,  the  expression  for  the  kinetic  energy  of  the  kin  mode  excited  by  a propagating  fault  is 
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and  for  a modulated  line  source: 


Ek=I  Z N0(W):Skm(V./2)| 

m=-i  ^ -4xm 


(III-7b) 


It  is  obvious  that  when  the  fault  length  approaches  zero,  the  expression  above  becomes 
identical  to  Eq.  ( III  -la);  also,  when  $ — 0,  $rQ/v  = T,  this  expression  is  identical  to  Eq.(III-lb). 
Similarly,  application  of  Eq.  (III-7b)  will  lead  to  Eq.  (Ill - 1 c)  under  these  conditions.  Nc^te  that 

, +0O 

(III—  8) 


27  f I CR(w)  | dw  = 1 
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thus,  the  kinetic  energy  of  a mode  is  independent  of  attenuation. 

The  results  for  a modulated  line  source  are  presented  in  Table  III- 2 and  in  Figs.IIf-15 

-2 

to  III-17.  Our  choice  of  the  model  has  been  dictated  by  its  asymptote  go  amplitude  behavior 
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at  high  frequencies,  cf.  Eq.  (III-3b),  in  agreement  with  Aki.  These  are  all  for  the  I960  Chilean 

1 7 

earthquake  using  Plafker  and  Savage's  fault-plane  solution  with  the  line  source  at  79.4-km 
depth.  Table  III-2  lists  the  cumulative  energy  in  the  toroidal  modes  for  a variety  of  fault-length 
and  rupture- velocity  combinations.  Thus,  one  can  see  here  the  effects  of  three  experiments: 

(1)  changing  the  fault  length  with  a constant  rupture  velocity  (3.5  km/sec);  (2)  changing  the  rup- 
ture velocity  with  a constant  fault  length  (800  km);  and  (3)  changing  the  fault  length  and  rupture 
velocity  but  keeping  the  ratio  constant  (180  sec).  The  most  rapid  variation  in  the  energy  values 
given  in  Table  III-2  occurs  when  the  rupture  velocity  is  comparatively  high:  3.5  to  4.0  km/sec. 

Figures  III- 1 5 to  III  - 1 7 show  typical  toroidal-energy  spectral  distributions  from  the  three 
above-mentioned  experiments.  In  particular.  Fig.  Ill- 1 5 results  from  a low  rupture  velocity 
(2.0  km/sec).  There  the  spectral  dropoff  rate  is  similar  to  the  half- sine -wave  point  sburce, 
but  the  scallops  in  the  power  envelope  are  not  as  prominent.  Figure  III- 1 6 shows  the  results 
from  a high  rupture  velocity  (4.0  km/sec).  These  are  similar  to  those  for  a step-function  point 
source.  This  is  explained  by  the  fact  that  a moving  source  radiates  energy  most  efficiently  in 
its  direction  of  travel,  and  is  spectrally  enhanced  at  high  frequencies  by  the  Doppler  effect. 

Seen  from  the  front,  where  the  radiated  energy  is  concentrated,  it  would  appear  as  a step- 
function  point  source.  Finally,  Fig.  Ill- 1 7 shows  the  results  from  intermediate  values  of  fault 
length  and  rupture  velocity.  Here,  the  energy  distribution  is  smooth  with  a peak  near  the 
160-sec  period.  The  distribution  also  drops  off  cleanly  above  this  peak. 

Results  for  the  line-source  model  with  finite  rupture  velocity  do  not  substantially  £lter  our 
previous  conclusions.  In  particular,  the  cumulative  toroidal  energy  which  we  computejd  for  the 

I960  Chilean  earthquake  using  the  fault-length  and  rupture-velocity  values  estimated  by  Kanamori 
18  "23 

and  Cipar  (800  and  3.5  km/sec)  is  9.1  X 10  ergs.  This  corresponds  to  a ha  If -sine  ■►wave 
point-source  duration  of  140  sec,  which  is  very  close  to  that  value  obtained  when  the  sburce  du- 
ration was  adjusted  to  give  the  reported  M value.  We  can  therefore  conclude  that,  for  the  I960 
Chilean  earthquake,  both  point-  and  line- source  models  give  energies  in  good  agreement  with 
the  Gutenberg-Richter  relation. 

The  principal  difference  between  the  point-  and  line-source  results  is  in  the  shape  of  the 
power  spectrum  envelope.  Going  to  a more  realistic  finite  source  with  a propagating  rupture 
gives  a smoother  power  spectrum.  However,  increasing  the  rupture  velocity  to  4.0  km/sec  from 
3.5  km/sec  produces  an  unacceptable  envelope  shape  which,  in  our  frequency  band,  resembles 
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TABLE  II 1-2 

TOTAL  TOROIDAL  ENERGY  FOR  THE  1960  CHILEAN  EARTHQUAKE1- 
AS  A FUNCTION  OF  FAULT  LENGTH  AND  RUPTURE  VELOCITY 

Toroidal  Energy 

Fault  Length 

Rupture  Velocity 

Ratio 

(erg) 

(km) 

(krry/sec) 

(sec) 

24 

8.96 X 10 

200.0 

3.5 

57.1 

24 

3.98 X 10 

400.0 

3.5 

114.3 

24 

1 .85  X 10^4 

600.0 

3.5 

171.4 

9.10X  1023 

800.0 

3.5 

228.6 

4.84  X 1023 

1000.0 

3.5 

285.7 

1 .09  X 1023 

800.0 

2.0 

400.0 

1 .09  X 1023 

800.0 

2.5 

320.0 

3.07  X 1023 

800.0 

3.0 

266.7 

9. 10  X 1023 

800.0 

3.5 

228.6 

24 

2.51  X 10 

800.0 

4.0 

200.0 

3.01 X 1023 

7.2 

0.04 

180.0 

2.70 X 1023 

36.0 

0.2 

180.0 

2.74  X 1023 

72.0 

0.4 

180.0 

4.48 X 10 

360.0 

2.0 

180.0 

24 

2.96  X 10 

720.0 

4.0 

180.0 

t Source  at  79.4-km  depth,  Plafker  and  Savage's^  fault-plane  parameters. 
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that  due  to  a step-function  point  source.  Since  the  shear  velocity  of  the  1066 A model  at  80 -km 
depth  is  approximately  4.5  km/sec,  we  conclude  that,  for  this  event,  the  radiated  energy  is 
quite  sensitive  to  rupture  propagation  velocity  approaching  the  source  region  shear  veloicity. 

D.  W.  McCowan 
A.  M.  Dziewonski 

G.  RAYLEIGH-WAVE  EXPLOSIVE  MOMENTS  AND  SEISMIC  ENERGIES 
FOR  SIX  U.S.  EXPLOSIONS 

1 9 

Rayleigh-wave  point-source  theory  provides  a means  of  calculating  the  radiated  Rjayleigh- 

e case 
vill  be: 


wave  seismic  energy  if  the  moment  tensor  and  time  history  of  a source  are  known.  In  tl 
of  a step  function  of  explosive  moment,  the  vertical  component  Rayleigh-wave  spectrum 

a/  % MkAk  y , v f~\  r , y , „ , dw(w.k.d)1 

A(r,z,u;)  j w(o>,k,z)  [-ku(u;,k,d)  + — * ■ ] 

ZTUjJ 


(III-9) 


k is 


In  this  expression,  A is  the  zero-to-peak  amplitude,  M is  the  scalar  explosive  moment, 

horizontal  wavenumber,  o>  is  angular  frequency,  d is  the  source  depth,  and  u and  w ard  the 

1 9 

Rayleigh-wave  eigenfunctions  normalized  as  described  previously.  In  the  interests  of  sim- 
plicity, the  amplitude  effects  of  attenuation  and  dispersion  have  been  neglected.  Figures  till  - 1 8 

23 

to  III-20  show  this  function  plotted  for  the  case  of  a 10  dyn-cm  explosive  moment  soured  at  a 
distance  of  100  km.  Figures  III- 1 8 and  III-19  show  the  effects  of  a change  in  source  depth, 
while  Fig.  Ill- 1 9 shows  the  results  of  Fig.  Ill- 1 8 seen  through  an  LRSM  (Long  Range  Seismic 
Measurements)  LP  recording  system. 

Values  read  from  these  graphs  or  similar  versions  computed  for  the  required  source  depths 
can  be  used  to  estimate  the  explosive  moment  of  an  underground  explosion  from  amplitude 
period  observations.  The  formula  is: 


and 


M = 


_ A(obs)  n/D 
A(graph)  10 


X 1023  dyn-cm 


(III- 10) 


where  D is  the  source  distance  in  kilometers.  Combining  Eq.  (Ill- 1 0)  with  the  Prague  formula 
for  Mg,  and  using  the  value  in  Fig.  III-18  at  20 -sec  period,  implies  the  following  surface-wave 
magnitude-moment  relation: 


Mg  = log  (M  x 10 


-22 


) + l.l6  logA  + 0.80 


(IIfl-11) 


where  A is  the  epicentral  distance  in  degrees,  and  M is  again  the  explosive  moment. 

2 0 

The  arguments  used  previously  to  derive  an  energy-moment  tensor  relation  appropriate 
to  free  oscillation  seismology  can  be  easily  generalized  to  Rayleigh  waves.  The  resulting  ex- 
pression is 


E = i Z J0  kdk  E 


m=-2 


M:S“ 

2lrwj2lR^wj’k* 


(III-12) 


which  is  valid  for  a step- source  time  function.  Here,  M is  the  moment  tensor  of  the  source 

R + 

Sm  is  the  set  of  Rayleigh-wave  strains,  1^  is  the  normalization  integral,  and  j is  the  mode 
index.  The  Rayleigh-wave  strains  can  either  be  calculated  by  using  a surface-wave-dispersion 
program  or,  for  a homogeneous  halfspace,  be  expressed  analytically.  In  the  latter  case21  tl)e 
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TABLE  111-3 

PUBLISHED  INFORMATION  ON  NTS  EXPLOSION  AND  COMPUTED  RESULTS 

" 'Shot(1) 

Date(1) 

Location^ 

Depth^ 

(km) 

r i 0) 
Geology 

Yield*” 

(kl) 

Explosive  Moment 

(dyn-cm)x  10^ 

(2) 

Step-Function' 
Energy  (kT) 

Scale(3) 

Length 

(km) 

RiseW 

Time 

(sec) 

Half  Beak 

6/30/66 

Pahute  Mesa 

0.88 

Rhyolite 

300 

49  ± 32 

400.0 

4.6 

1.4 

Greeley 

] 2/20/66 

Pahute  Mesa 

1.23 

Tuff 

825 

165  db 

190 

1800.0 

8.0 

2.4 

Scotch 

5/23/67 

Pahute  Mesa 

1.00 

Tuff 

150 

33  ± 

29 

160.0 

4.4 

1.3 

Knickerbocker 

5/26/67 

Pahute  Mesa 

0.63 

Rhyolite 

71 

20  ± 

16 

200.0 

4.7 

1.4 

Gasbuggy 

iyio/67 

Dulce,  NM 

1.29 

Shale 

29 

3.6  ± 

2.2 

0.73 

0.77 

0.23 

Rulison 

9/10/69 

Rulison,  CO 

2.57 

Sandstone 

40 

3.3  ± 

1.7 

0.08 

- 

- 

^Taken  from  Ref. 23. 

<2,1  kl  = 4 X 
<3)-ct. 

in19 

10  ergs. 

^Assuming  c = 

= 3.3  krr/sec. 

resulting  energy  density  consists  of  three  terms,  of  the  form: 


2 -Atcd 
w e 


(III- 13) 


Each  one  of  these  decays  at  both  a;  extremes.  The  energy  for  a homogeneous  Poisson  halfspace 
can  be  integrated  to  give  an  exact  answer  for  the  radiated  energy: 


r,  (0.1  96)  M 

E 13 


l±d 


III-14) 


22  24  25 

However,  previous  free-oscillation  results  ' ' show  that  a step  function  in  applied  moment 

produces  an  infinite  radiated  energy.  Equation  (III-l 4)  shows  that,  for  fundamental-mode  Ray- 
leigh waves,  this  only  happens  at  zero  source  depth.  For  nonzero  source  depths,  the  exponen- 
tial in  Eq.  (Ill- 1 3)  makes  the  energy  finite.  In  other  words,  very-high-frequency  Rayleigh  waves 
propagate  between  the  free  surface  and  the  source  and  are  therefore  poorly  excited.  Instead  of 
relying  on  this  source-depth  effect  to  limit  the  high-frequency  energy  spectrum,  it  is  prucfent  to 
include  the  possibility  of  a ramp-source  time  function.  The  corresponding  expressions  to 
Eqs.  (in-12)  and  (III- 1 4)  are: 


+2 

=?  i r 


m=-2 


kdk  Z 

3 


|M:S*  («  k,d)|' 
— J 

Znw.  IR(w.,k) 


sin  (ukT/2) 


(sin  1 


TfJZ 


(I1I-15) 


and 


E = 


(0.196)M2 

2.2909 

3.0475  ! 1.0343 

Hd3 

1 4.8698  + 1.6950x2 

1.9103  + 1.2408x2  0.4867  + 0.7866x2. 

x = cT/d 


(III-16) 


where  T is  the  ramp  time.  The  bracketed  quantity,  denoted  by  2(x),  is  plotted  in  Fig.  III-31. 

As  can  be  seen,  increasing  the  ramp  time  serves  to  substantially  decrease  the  radiated  energy. 

The  results  of  estimating  explosive  moments  and  computing  radiated  energies  from  Eqs.  (Ill -9 ) 

and  (III— 1 6 ),  respectively,  are  given  in  Table  III- 3 . These  are  for  a selection  of  six  U.S.  explo- 

25-3  0 

sions  described  as  part  of  the  LRSM  program.  In  applying  Eq.  (Ill- 9).  a structure  appropri- 

ate to  the  NTS  was  used?1  The  explosive  moments  and  their  standard  deviations  were  computed 
from  all  the  nonquestionable  Rayleigh-wave  data  reported  for  these  six  explosions.  No  attempts 
were  made  to  remove  anomalously  high  or  low  values  (which  were  clearly  evident),  nor  wer^ 
any  attempts  made  to  reduce  the  impact  on  the  average  of  very  distant  stations.  A value  of 
(1  = 4 X 10d 
analyses. 

The  energies  in  Table  M-3  are  on  the  order  of,  or  slightly  higher  than,  the  announced  yields 

for  NTS  explosions.  The  two  gas  stimulation  explosions,  Gasbuggy  and  Rulison,  produce  energy 

estimates  much  lower  than  the  announced  yields.  These,  of  course,  are  all  for  the  step-function 

source.  In  an  ad  hoc  attempt  to  deduce  the  correct  source  ramp  time,  scale  lengths  (=cT)  were 

chosen  which  reduced  the  calculated  step-function  energy  to  1 percent  of  the  announced  yield. 

32 

This  was  thought  to  be  a representative  number  based  on  previous  work.  The  resulting  ramp 
times  are  all  near  1.5  sec  for  the  NTS  events,  except  that  for  the  Greeley  explosion.  Unfortu^ 
nately,  its  explosive  moment  is  the  most  poorly  determined  example  in  the  group.  Again,  the 


dyn/cm2,  deemed  characteristic  of  sedimentary  crustal  rocks,  was  used  in  all  ihe 
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TABLE  111-4 


PARAMETERS  FOR  pP  AND  SLAPDOWN  PHASES  DETERMINED  BY  (a)  BAKUN  AND  JOHNSON34 
AND  (b)  KING  ET  AL.?5  CORNER  FREQUENCIES  WERE  CALCULATED  FROM  KING'S  DATA. 

mb  AND  Ms  ARE  NEIS  MAGNITUDES. 


Amplitude  Ratio 

Time  Delay  (sec) 

Comer  Frequency 
(Hz) 

Yield 

(kT) 

pp/p 

Slapdown/P 

pp/p 

Slapdown/P 

mb 

M 

s 

Longshot  (b) 

-0.3 

+0.3 

0.5 

0.9 

1.62 

85 

6.1 

3.9 

(a) 

- 0.7 

+0.7 

0.8 

1.4 



_ 





Milrow 

(b) 

-0.5 

+0.4 

0.8 

1.4 

1 .21 

-1000 

6.5 

4.9 

(a) 

—0.6 

+0.8 

1.2 

1.9 









Cannikin 

(b) 

-0.4 

+0.4 

1.1 

1.9 

0.971 

5000 

6.8 

5.7 

gas  stimulation  explosions  stand  out  as  having  substantially  lower  ramp  times.  In  fact  the 
computed  step-function  energy  for  Rulison  is  already  below  1 percent  of  the  announced  ydeld. 

An  obvious  conclusion  of  this  work  is  that  the  source  time  function  for  explosions  remains 
virtually  unknown.  Observed  Rayleigh-wave  spectra  corrected  for  the  recording  systerji  re-’ 
sponse30  generally  show  that  the  energy  peaks  above  l/lO  Hz.  Since  this  is  in  the  rejection 
band  of  the  anti-microseism  filter,  there  isn’t  much  chance  of  recovering  the  source  tirhe  func- 
tion from  presently  available  data.  What  is  needed  is  a representative  sample  of  broadband 
Rayleigh-wave  seismograms  recorded  on  reasonably  near-in  stations.  One  good  possibility 
here  is  to  modify  the  Albuquerque  SRO  by  reducing  the  gain  and  removing  the  anti-micrqseism 


filter.  It  could  then  be  used  to  observe  NTS  explosions  which,  by  chance,  propagate  to  Albu- 
querque across  a relatively  simple  and  homogeneous  structure:  the  Colorado  Plateau.  These 
recordings  could  then  be  analyzed  along  the  lines  of  Dziewonski  and  Gilbert's*^  work  on  (Jeep 
South  American  earthquakes.  The  result  of  such  endeavors  should  be  an  accurate  determination 
of  the  explosive  moment  time  histories  for  NTS  explosions,  which  may  also  lead  to  more  ac- 
curate estimates  of  the  radiated  seismic  energies. 

D.  W.  McCowan 

H.  YIELD  VERIFICATION  USING  LONG-PERIOD  P WAVES  FROM  EXPLOSIONS 

At  WWSSN  stations,  intermediate-  to  large-yield  nuclear  explosions  produce  LP  P whves 
which  are  extremely  simple  in  shape.  Although  these  waves  are  not  usually  seen  for  yieldjs 
much  less  than  100  kT,  their  simple  shape  and  generally  smaller  amplitude  scatter  at  large 
yield,  compared  with  SP  records,  allow  an  alternate  way  to  examine  explosion  source  functions 
and  seismic  scaling  with  yield. 

Figures  III-22(a),  (b),  and  (c)  show  plots  of  LP  and  SP  amplitudes  in  microns  for  three  large 
explosions  in  different  test  areas.  These  amplitudes  were  measured  peak-to-peak  on  available 
WWSSN  records,  and  corrected  only  for  instrument  gain  at  1 sec  for  SP  stations  and  at  15  sec 
for  LP  stations.  Except  for  the  Novaya  Zemlya  data,  in  Fig.  III-22(a),  the  LP  amplitudes  are 
mostly  greater  than  the  SP  amplitudes,  with  a slope  LP  vs  SP  less  than  1.  For  both  the  Ncivaya 
Zemlya  and  Cannikin  tests,  the  LP  amplitudes  have  an  order-of-magnitude  less  scatter  than  the 
SP  amplitudes. 

It  is  interesting  that  the  LP  amplitudes  for  Cannikin  are,  on  the  average,  about  half-a- 

magnitude  larger  than  the  Novaya  Zemlya  LP  amplitudes,  although  the  m^  values  assigned  b^y 

the  NEIS  were  the  same.  This  observation  and  the  larger  Cannikin  yield  suggest  that  the  Cah- 

nikin  m,  is  underestimated, 
b 

In  order  to  explain  the  pattern  of  SP  and  LP  amplitudes,  we  investigated  the  effect  of  ext 
plosion  source  functions  and  earth  attenuation  on  synthetic  SP  and  LP  records.  Because  our 
knowledge  of  site  parameters  for  Novaya  Zemlya  explosions  is  limited,  we  chose  the  three  UlS. 
Amchitka  tests  - Longshot,  Milrow,  and  Cannikin  - which  have  a wide  range  of  known  yields 
and  depths. 

Table  III-4  contains  source  parameters  for  pP/P  and  slapdown/P  phases  obtained  by  Bakun 
34  35 

and  Johnson  and  by  King  et  al.  from  SP  spectral  analysis.  - 

Although  the  delay  times  for  pP  and  slapdown  phases  relative  to  P agree  between  the  two 
studies,  the  relative  amplitudes  do  not.  In  particular,  pP  and  slapdown  phases  have  about 
equal  and  opposite  amplitudes  significantly  less  than  1. 
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Synthetic  LP  WWSSN  records  for  Cannikin  were  calculated  to  compare  with  actual  wave- 
forms. Included  in  the  records  are  Bakun  and  Johnson’s  parameters  for  the  relative  P,  pP,  and 
slapdown  amplitudes,  an  explosion  source  function  given  by  Blake^  with  a corner  frequency  of 
0.971  obtained  from  King's  results,  a causal  attenuation  operator,  and  the  WWSSN  LP  response. 

Figure  III-23(a)  shows  the  synthetic  waveforms  for  the  P phase  alone,  and  the  effects  of 
adding  pP  and  slapdown  phases.  Figure  III-24  displays  several  LP  recordings  of  Cannikin,  and 
they  are  surprisingly  impulsive.  Of  particular  interest  are  the  equal  positive  and  negative 
swings  in  the  first  motion,  and  the  lack  of  strong  later  phases  in  the  first  few  seconds.  This 
seems  to  imply  that  the  pP/P  amplitude  ratio  must  be  closer  to  one,  and  the  slapdown  ampli- 
tude nearer  to  zero  than  suggested  by  Bakun  and  Johnson's  study.  The  slapdown  phase  shown 
in  the  synthetic  record  is  not  visible  at  all  on  most  of  the  LP  explosion  records. 

In  Fig.  Ill— 2 3 (b) , synthetic  LP  records  for  Cannikin  are  displayed  assuming  a pP/P  ampli- 
tude ratio  of  -1.0  and  retaining  the  slapdown/P  amplitude  ratio  of  0.8.  The  synthetic  records 
which  contain  only  P and  pP  phases  adequately  explain  the  shape  and  pulse  width  of  the  observed 
data.  In  Fig.  III-23(c),  synthetic  SP  records  are  shown  for  Cannikin  using  the  same  source  pa- 
rameters. These  calculations  use  the  formulation  for  the  WWSSN  SP  response  obtained  by 
37 

Burdick  and  Mellman,  which  includes  the  effect  of  transducer  inductance  at  high  frequencies. 

Synthetic  SP  and  LP  records  were  computed  for  Cannikin,  Milrow,  and  Longshot  for  a 
variation  of  attenuation  parameter  t*  = (travel  time)/ (average  Q).  The  far-field  displacement 
spectrum  for  each  explosion  was  normalized  to  equal  the  yield  at  0 frequency,  and  given  the 
corner  frequency  listed  in  Table  III-4.  Instrument  responses  were  normalized  to  1.0  at  1 and 
15  sec  for  SP  and  LP  records,  respectively.  P amplitudes,  peak-to-trough,  were  measured 
to  obtain  relative  amplitudes  of  LP  and  SP  phases. 

These  synthetic  amplitudes  are  displayed  in  Fig.  III-2  5.  The  trend  of  SP  vs  LP  amplitudes 
is  very  similar  for  the  three  explosions  in  spite  of  variations  in  source  corner  frequency.  This 
trend  is  controlled  primarily  by  attenuation  variations  which  affect  the  SP  amplitudes  more 
severely  than  the  LP  amplitudes.  Thus,  attenuation  alone  may  explain  much  of  the  trend  of  the 
real  LP  and  SP  data  in  Figs.  III-22(a-c).  In  particular,  the  Novaya  Zemlya  data  show  a cross- 
over region  where  LP  amplitudes  become  less  than  SP  amplitudes.  From  the  synthetic  cases 
this  would  be  expected  along  paths  with  t*$,  0.5,  a reasonable  value  for  high  Q teleseismic 
paths. 

Another  feature  of  the  synthetic  data  is  that  attenuation  differences  may  cause  SP  ampli- 
tudes to  overlap  from  explosion-to-explosion,  whereas  the  LP  amplitudes  remain  separated. 
This  would  decrease  the  reliability  of  SP  magnitudes  based  on  uneven  station  coverage. 

Taking  amplitude  ratios  of  the  synthetic  data  for  pairs  of  shots,  we  obtain  the  curves  shown 
in  Fig.  III-26(a)  as  a function  of  t*.  As  attenuation  increases,  both  LP  and  SP  amplitude  ratios 
converge  to  the  ratio  of  yields  used  in  the  calculations. 

Figure  III-26(b)  shows  amplitude  ratios  calculated  from  real  data  of  various  periods.  On 
LP  records  we  measured  32  Cannikin  LP  amplitudes,  and  8 Milrow  amplitudes.  Of  these  mea- 
surements, Cannikin  and  Milrow  were  measured  at  7 common  stations  including  COL  which  also 
recorded  Longshot.  Amplitude  ratios  at  common  stations  are: 

Milrow/Cannikin  0.26  ± 0.07 

Longshot/Milrow  0.089 

Longshot/Cannikin  0.023 
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These  ratios  are  plotted  in  Fig.  III-26(b)  as  crosses.  In  spite  of  the  small  number  of  measure- 
ments, these  ratios  agree  well  with  the  yield  ratios  and  with  ratios  expected  from  NEIS  Mg 
differences  which  are  measured  at  periods  from  10  to  20  sec.  Amplitude  measurements  pre- 
dicted from  NEIS  m^  differences  clearly  do  not  agree  with  the  other  data  and  must  be  considered 
unreliable.  Specifically,  for  both  Milrow  and  Cannikin  are  probably  underestimated.  This 
is  suggested  by  SP  measurements  at  seven  common  LRSM  stations  by  von  Seggern  and  Blandford 
They  computed  the  amplitude  ratios  displayed  as  solid  points  in  Fig.  III-26 (b) . These  radios  are 
much  closer  to  the  yield  ratios  than  are  the  ratios  based  on  NEIS  m^. 

This  preliminary  analysis  of  LP  data  from  explosions  is  promising.  We  intend  to  examine 
data  from  all  test  areas  and  concentrate  on  the  yield  range  of  100  to  1000  kT. 
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ANOMALOUS  mb  VS  YIELD  FOR  LONGSHOT,  MILROW,  AND  CANNIKIN 

Using  LP  and  SP  data,  we  are  investigating  the  difficulties  of  estimating  consistent  yields 
for  large  explosions.  In  Sec.  H above,  we  indicated  that  relative  LP  amplitudes  for  body  jvaves 
from  Longshot,  Milrow,  and  Cannikin  could  be  used  to  estimate  the  relative  yields  of  the  three 
explosions.  These  amplitudes  were  not  divided  by  their  apparent  period,  nor  were  they  corrected 
for  instrument  response. 

Figure  III-27  shows  log  amplitude  as  a function  of  yield  for  Milrow  and  Cannikin  data,  nor- 
malized by  the  Longshot  data.  Different  types  of  actual  and  theoretical  amplitudes  are  shown, 

-2 

both  LP  and  SP.  Using  the  scaled,  theoretical  to  source  spectra  of  the  previous  section, 
calculated  amplitudes  of  Milrow  and  Cannikin  relative  to  Longshot  at  0.05,  1.0,  and  2.0 
quencies.  In  the  broadband  from  0.05  to  1.0  Hz,  the  theoretical  ratios  vary  insignificantly1  and 


we 

fre- 


M 


s' 

3 sec. 


agree  quite  well  in  Fig.  HI-27  with  amplitude  variations  expected  from  differences  in  NEIS 
measured  from  15  to  20  sec,  and  LP  body  wave  measurements  at  periods  as  small  as  2 tc 
Therefore,  a consistent  trend  of  amplitude  proportional  to  yield  1.0  fits  the  data  from  periods 
£,2  to  20  sec. 

At  shorter  periods  ~ 1 sec  the  amplitudes  are  proportional  to  a smaller  power  of  yield.  For 

example,  at  2 Hz  the  theoretical  amplitudes  vary  as  yield  ~0.75  because  the  corner  frequencies 

assumed  in  the  source  functions  for  Milrow  and  Cannikin  are  both  <2.0  Hz.  Von  Seggern  and 
38 

Blandford  measured  the  amplitudes  of  the  first  quarter-cycle  of  the  SP  records  at  seven 
LRSM  stations.  They  obtained  the  following  ratios: 

Milrow/ Longshot  6.68  ± 0.86 

Cannikin/ Longshot  15.76  ± 2.32 

Logarithms  of  these  ratios  are  shown  in  Fig.  III-27  as  black  circles,  and  follow  the  trend  of  the 
theoretical  curve  at  2 Hz.  Although  the  first  motion  amplitudes  contain  frequencies  higher  than 
1 Hz,  it  is  uncertain  whether  the  spectra  peak  nearer  to  2 and  1 Hz.  If  not,  the  data  of  von  Seggern 
and  Blandford  may  indicate  a faster  falloff  than  to  ^ in  the  source  spectra  beyond  the  corned 
frequencies. 

Expected  log  amplitudes  based  on  NEIS  and  ISC  m^,  assuming  similar  periods  for  each  P 
wave,  show  a trend  well  below  von  Seggern  and  Blandford' s results  and  must  be  considered  sus- 
pect. Cannikin  and  Milrow  m^  values  are  most  certainly  underestimated,  due  to  a variety  c|f 
problems  in  network  detection,  station  bias,  and  inconsistent  period  measurements.  These  are 
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discussed  in  great  detail  by  Basham  and  Horner.  They  used  the  Canadian  seismic  network  to 
estimate  LP  and  SP  magnitudes  which  do  not  involve  dividing  by  pulse  period.  By  including 
station  biases  for  high-quality  stations  on  the  Canadian  shield,  Basham  and  Horner  calculated 
SP  magnitudes  which  give  the  log  amplitude  ratios  shown  in  Fig.  III-27  by  black  triangles.  These 
values  agree  well  with  von  Seggern  and  Blandford’s  amplitude  ratios  and  verify  the  poor  quality 
of  NEIS  m,  values. 

k 40 

Basham  and  Horner  used  an  explosion  source  model  by  Haskell  to  predict  the  LP  and  SP 

body  phase  amplitudes  and  the  relative  Mg  values  for  the  three  Aleutian  explosions.  As  in  our 
study,  they  were  able  to  match  Mg  variations  based  on  the  0-Hz  spectral  level,  which  scales  as 
yield  1.0.  They  could  not,  however,  predict  either  the  LP  or  SP  amplitudes  using  the  Haskell 
source  model.  Their  discrepancy  for  SP  amplitudes  is  similar  to  that  shown  in  Fig.  III-27,  be- 
tween the  theoretical  curve  for  1.0  Hz  and  the  black  triangles  and  circles.  At  long  periods,  the 
scaled  Haskell  source  functions  predict  a Cannikin/Longshot  P-wave  ratio  about  0.4  magnitude 
units  larger  than  the  prediction  in  Fig.  III-27. 

One  difficulty  with  the  Haskell  explosion  model  is  that  it  predicts  a far-field  spectrum  which 
-4 

decays  as  co  above  the  corner  frequency.  This  causes  large  explosions  to  have  less  high- 

3 8 -2 

frequency  energy  than  small  explosions.  Most  evidence,  however,  suggests  that  cj  source 
models  fit  the  explosion  data  better  and  are  physically  more  plausible  than  Haskell’s  model. 

It  appears  then  that  m^  as  currently  determined  by  the  NEIS  and  ISC  is  particularly  ques- 
tionable for  large  explosions,  and  that  more  consistent  yield  determinations  can  be  made  by  the 
combined  use  of  Mg  and  LP  P waves.  It  must  be  emphasized  that  if  absolute  rather  than  rela- 
tive yields  are  to  be  obtained,  a path -dependent  magnitude  M as  proposed  by  Marshall  and 

41  s 

Basham,  must  be  used  - particularly  for  explosions  in  different  tectonic  regions. 

C.  W.  Frasier 
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Fig.  Ill— 1 . Amplitudes  vs  azimuth.  All  signals  were  recorded  by  LP  instruments 
of  WWSSN  in  distance  range  30°  to  90°  and  reduced  peak-to-peak  amplitudes  in  mi- 
crons normalized  to  an  arc  distance  of  60°.  Corrections  for  spherical  spreading 
were  computed  from  earth  model  CIT20  (Ref.  5).  Events  used  are  (a)  Cannikin  ex- 
plosion, and  (b)  Novaya  Zemlya  event  (14  October  1970). 
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Fig.  III-2.  Relocated  epicenters  of  Kurile  aftershock  sequence 
of  21  January  - 10  February  1976. 
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Fig.  Ill- 3.  LP  vertical  Rayleigh -wave  trains  recorded  at  Mashad 
SRO  from  Kurile  sequence,  covering  3-l/2  orders  of  magnitude. 
All  seismograms  start  at  2070  sec  after  origin  time,  and  are  of 
640-sec  duration. 
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Fig.  HI-4.  Spectra  of  seismograms 
(a)  through  (d)  from  Fig.  III-3. 
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together  with  ratios  predicted  by  a; 
law. 
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Fig.  III-7.  Final  averaged  WWSSN  SP  calibration  pulse  recorded 
at  Weston  Seismic  Observatory. 
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Fig.  III-8.  System  amplitude  response  computed  from  pulse 
in  Fig.  III-7.  F 


Fig.  III-9.  System  phase  response  computed  from  pulse 
in  Fig.  Ill -7. 
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FREQUENCY 


Fig.  Ill-lO(a-c).  Observed  variation  of  frequency  with  log  A/T  for  three  stations,  using 
ISC  Bulletin  data.  Solid  line  is  an  inferred  seismicity  relation,  and  dashed  lines  are  in- 
terpretations of  station  detection  capability.  Gd  is  50 -percent  detection  threshold;  *Gg  is 
50 -percent  saturation  threshold. 
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Fig.  Ill-ll(a-c).  Observed  variation  of  frequency  with  station  magnitude  mb,  for  events  in 
distance  range  30°  to  90°,  using  ISC  Bulletin  data.  For  other  symbols,  see  Fig.  III-10. 
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Fig.  Ill- 1 Z.  Comparison  of  magnitude  distri- 
bution of  events  input  into  simulated  network 
program,  with  magnitude  distribution  of  events 
detected.  Gs  and  Gd  are  50 -percent  detection 
and  saturation  thresholds  of  network. 


Fig.  III-13.  Comparison  of  true  magni- 
tudes of  events  detected  by  simulation 
program  with  network  magnitudes,  com- 
puted as  average  station  magnitude  for 
those  stations  detecting  each  event. 
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Fig.  Ill -14.  Average  values  of  difference  between  network  average  magnitude 
and  true  magnitude,  plotted  as  a function  of  true  magnitude. 
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Fig.  m-15.  Spectral  distribution  of  radiated 
energy  for  I960  Chilean  earthquake  for  an 
800 -km  fault  with  2.0-km/sec  rupture  ve- 
locity. Line  source  at  79.4-km  depth  using 
Plafker  and  Savage's17  fault-plane  solution. 


Fig.  Ill- 1 6 . Spectral  distribution  of  radiated 
energy  for  I960  Chilean  earthquake  for  a 
720 -km  fault  with  4.0-km/sec  rupture  ve- 
locity. Line  source  at  79.4-km  depth  using 
Plafker  and  Savage's17  fault-plane  solution. 
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Fig.  Ill- 1 7 . Spectral  distribution  of  radi- 
ated energy  for  I960  Chilean  earthquake 
for  an  800 -km  fault  with  3.5-km/sec  rup- 
ture velocity.  Line  source  at  79. 4 -km 
depth  using  Plafker  and  Savage’s1  ' fault- 
plane  solution. 
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Fig.  HI-18.  Amplitude  vs  period  for  a 
23 

10  dyn-cm  explosive  moment  source 
1 km  deep  at  100-km  distance.  (Instru- 
ment response  not  included.) 
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Fig.  Ill— 1 9-  Amplitude  vs  period  for  a 
23 

10  dyn-cm  explosive  moment  source 
3 km  deep  at  100-km  distance.  (Instru- 
ment response  not  included.) 


Fig.  III-20.  Amplitude  vs  period  for  a 

23 

10  dyn-cm  explosive  moment  source 
1 km  deep  at  100-km  distance  as  seen 
through  LRSM  LP  recording  system. 


Fig.  Ill -21.  2 factor  for  ramp-source 

time  function  in  a homogeneous  Poisson 
halfspace. 


72 


LP  AMPLITUDE  (/i) 


Fig.  m-22.  LP  and  SP  P amplitudes  measured  at  WWSSN  stations  for  three  large 
explosions  at  (a)  Novaya  Zemlya,  (b)  Amchitka  Island,  and  (c)  NTS. 


73 


Fig.  III-23.  Synthetic  P waves  from  Cannikin.  Includes  explosion  source  spectrum,  earth 
attenuation  pP  and  slapdown  phases,  and  WWSSN  response,  (a)  LP  synthetic  record  with 
amplitude  ratios  pP/P  = -0.6,  slapdown/P  = 0.8,  and  t*  = 0.5.  (b)  LP  synthetic  records 

with  pP/P  = -1.0,  slapdown/ P = 0.8,  and  three  attenuation  parameters.  (c)  Same  as  (b), 
except  recorded  on  SP  seismograph. 
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Fig.  III-24.  Examples  of  WWSSN  LP  recordings  of  P waves  from  Cannikin. 
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Fig.  III-25.  Relative  LP  and  SP  amplitudes  predicted  for  Cannikin,  Milrow, 
and  Longshot.  Source  spectra  have  corner  frequencies  given  in  Table  III-4, 
and  have  O-Hz  levels  scaled  by  yield  in  kilotons.  Dashed  lines  join  ampli- 
tude points  of  equal  attenuation  parameter  t*. 
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(o)  (b) 

Fig.  III-26.  (a)  Theoretical  amplitude  ratios  for  pairs  of  explosions  calculated 

from  Fig.  ni-25,  as  a function  of  t*.  (b)  Measured  LP  amplitude  ratios  (+),  SP 

ratios  determined  by  von  Seggern  and  Blandford38  (•),  yield  ratios  (*),  ampli- 
tude ratios  inferred  from  NEIS  Ms  (A)  and  (O). 
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Fig.  III-27.  Log  amplitude  as  a function  of  yield  for  Milrow  and  Cannikin 
seismic  phases  relative  to  Longshot.  Amplitudes  are  based  on  several 
SP  body-wave  magnitudes,  surface-wave  magnitudes,  and  some  LP  and 
SP  amplitude  measurements.  Curves  show  low  amplitude  variations  at 
0.05,  1.0,  and  2.0  Hz  predicted  by  an  explosion  source  model  de- 

scribed in  text. 
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rv.  MISCELLANEOUS  STUDIES 


A.  MULTIPLE-CHANNEL  MAXIMUM- ENTROPY  SPECTRA 

Multivariate  spectral  analysis  provides  a method  by  which  correlation  information  existing 
between  several  time  series  can  be  uncovered  and  evaluated.  For  seismic  data,  the  sha;>e  and 
amplitude  of  the  spectrum  of  each  seismogram  are  themselves  useful  information,  in  particular 
when  the  form  of  the  spectrum  is  predictable.  In  the  bands  where  the  various  spectra  are  strongly 
correlated,  multivariate  spectral  analysis  can  be  applied  to  deduce  the  relative  phase  and  group 
relationships.  The  success  of  such  spectral  analyses  in  this  respect  has  been  demonstrated  by 
the  use  of  the  maximum -likelihood  method  in  wavenumber  analysis  of  array  data. 

Conventional  methods  for  estimating  cross  spectral  density  are  based  on  generalizations  of 

ences. 

In  the  Blackman  and  Tukey  approach,  the  matrical  autocorrelation  function  is  estimated  {lirecUy 
from  the  data.  This  function  is  scaled  by  a taper  to  window  out  the  larger  lags  that  are  inaccu- 
rate due  to  finite  sample  size.  The  resulting  matrix  sequence  is  Fourier  transformed  tc  pro- 
duce the  matrical  spectral  sequence.  Similarly,  in  the  periodogram  approach,  window  functions 
are  employed  to  smooth  out  the  effect  of  finite  sample  size. 

These  approaches  as  applied  to  single -channel  data  have  shortcomings  that  do  not  disappear 
in  their  multivariate  extensions.  Specifically,  the  effect  of  windowing  is  to  smooth  the  p£wer 
spectrum  and  so  reduce  the  resolution.  Recently,  new  spectral  estimation  techniques  th^it  do 

not  make  assumptions  about  the  data  outside  the  observation  interval  have  been  shown  to  be  supe- 

4-7 

rior  in  terms  of  resolution  on  single -channel  data.  The  basic  method,  introduced  to  geophys- 

g 

ics  by  Burg,  is  based  on  maximizing  the  information  entropy  of  a time  series,  and  so  is  referred 

to  as  the  Maximum  Entropy  Method  (MEM).  Equivalently,  the  technique  can  be  thought  df  as 

9 

fitting  an  autoregressive  model  to  the  time  series,  and  is  referred  to  as  autoregressive  spectral 
analysis  in  the  engineering  literature. 

There  are  essentially  two  algorithms  for  computing  the  autoregressive  coefficients.  The 

10 

first  is  by  onestep  forward  prediction  error.  The  second,  developed  by  Burg,  combines  for- 
ward and  backward  prediction  to  insure  positive  definite  autocorrelation  functions.  In  thjis  re- 
port, we  investigate  the  extension  of  the  single -channel  Burg  technique  to  multichannel  data. 

Various  aspects  of  the  Multivariate  Burg  Spectral  Estimator  were  investigated.  Specifically, 
the  accuracy  and  resolution  of  estimated  spectra,  the  determination  of  optimal  autoregressive 
order,  and  the  effect  of  additive  white  noise  on  the  correctness  of  an  autoregressive  mo4el  were 

evaluated.  Although  theoretical  study  of  the  properties  of  an  autoregressive  spectral  estimator 
1 1 

has  been  attemped,  a complete  and  satisfactory  theory  has  not  been  developed.  The  investiga- 
tion undertaken  here  was  empirical.  Owing  to  the  large  variety  of  data  encountered  in  geophys- 
ics, no  statement  based  on  empirical  investigation  can  be  conclusive.  However,  studies)  based 
on  idealized  data  types  can  be  very  beneficial.  Other  than  illuminating  various  characteristics 
of  the  Multivariate  Burg  Estimator  under  a we  11 -controlled  and  adjustable  experimental  environ- 
ment, these  studies  provide  insights  in  understanding  more  complicated  spectra. 

In  general,  a spectrum  can  be  decomposed  into  broadband  and  narrowband  components.  Al- 
though the  boundary  between  these  is  not  well-defined,  the  flat  spectrum  for  a white  process  and 

the  line  spectrum  for  multiple  sinusoids  have  been  tested  to  examine  the  extreme  forms*  The 

12 

algorithm  for  computing  the  Burg  Maximum  Entropy  Spectra  was  taken  from  Nuttal. 
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First,  two-channel  100 -point  broadband  series  contaminated  with  1 -percent  white  noise 
and  with  a constant  time  shift  between  the  channels  were  analyzed  for  cross -spectral  amplitude 
coherence  and  phase.  Model  orders  equal  to  the  delay  were  used,  and  the  time  shift  was  kept 
small  compared  with  the  data  length.  It  was  observed  that,  in  general,  the  maximum  and  min- 
imum cross  powers  were  within  3 dB  of  the  true  value.  Values  of  the  estimated  coherence  at 
all  frequencies  were  well  above  0.97,  with  maxima  close  to  0.99,  which  is  the  true  coherence 
for  a signal-to-noise  ratio  (SNR)  of  100.  Phase  estimates  were  strikingly  linear,  varying  by 
less  than  3 percent  of  the  true  value  over  the  entire  spectrum. 

As  for  accuracy,  the  superiority  of  this  method  over  the  windowing  method  is  illustrated  in 
Figs.  IV-1  and  IV-2  where  coherence  and  phase  estimated  by  Burg's  method  are  compared  with 
those  computed  by  taking  the  Fourier  transform  of  a windowed  estimate  of  the  crosscorrelation 
sequence.  Two  coherence  and  phase  spectra  are  shown,  one  with  a 6 -sample  Bartlett  window, 
and  the  other  with  a 100-sample  Bartlett  window.  Even  the  100-sample  Bartlett  estimates  can- 
not reach  the  same  accuracy  achieved  by  the  Burg  estimates. 

The  same  procedure  was  also  applied  on  data  composed  of  two  sinusoids  with  constant  time 
shifts  between  channels  and  small  amounts  of  white  noise.  The  cross- spectral  amplitudes  showed 
the  two  spectral  peaks  to  be  well -resolved,  although  the  peak  values  varied  by  up  to  3 dB  from 
the  true  power  levels.  They  were,  however,  consistent  when  the  area  under  each  spectral  peak 
was  considered.  Coherences  at  the  frequencies  of  the  sinusoids  were  practically  equal  to  one. 
However,  they  were  not  as  well-resolved  as  the  cross  amplitude.  Phases  at  the  frequencies 
were  within  5 percent  of  the  true  value,  an  error  that  is  within  one  sampling  interval  and  so  quite 
acceptable.  As  was  previously  done  for  broadband  noise,  the  accuracy  of  the  Burg  method  is 
compared  with  windowing  methods  in  Figs.  IV-3  and  IV-4  for  an  input  series  of  10  cycles  of  a 
single  sinusoid  and  a 0.6  r shifted  version  in  the  second  channel.  The  Burg  estimate  of  order  4 
is  clearly  more  accurate  in  both  cross  amplitude  and  phase  than  the  Bartlett  estimators. 

Our  next  effort  involved  model-order  determination  for  multichannel  data.  The  Final  Pre- 

13  • 14  9 

diction  Error,  the  Information  Theoretic,  and  the  Autoregressive  Transfer  Function  criteria 

were  used  with  the  following  result:  the  order  numbers  predicted  by  each  method  generally 

o 

agree.  Exceptions  do  occur,  and  in  these  cases  use  of  Parzen1  s number  produces  the  best 
spectra.  For  processes  where  the  number  of  points  that  channels  are  shifted  is  greater  than 
the  order  of  the  individual  channels,  the  order  number  generated  is  the  shift  number.  The  re- 
sult is  better  cross  spectra  than  individual  spectra.  More  complex  individual  channels  result 
in  orders  larger  than  the  shifts. 

As  an  example  of  the  effect  of  order  number,  the  coherence  (Fig.  IV- 5)  and  phase  (Fig.  IV-6) 
for  the  two  shifted  white  processes  described  earlier  were  computed  by  both  the  Burg  and  win- 
dowing methods  as  functions  of  order  number.  Each  of  the  order-number  criteria  gave  the  shift 
(i.e.,  10)  as  the  order  number.  Burg's  estimate  rapidly  converged  to  the  true  spectra  at  and 
beyond  the  autoregressive  order  of  10.  Phase  estimates  from  the  other  method  are  comparable 
to  Burg's  estimates  only  when  autocorrelation  lags  of  20  and  longer  are  used.  Coherence 
estimates,  nevertheless,  seen  incapable  of  achieving  the  same  level  of  accuracy  as  Burg's 
estimates,  even  when  50  samples  of  the  autocorrelation  sequence  were  used. 

Our  final  effort  looked  at  the  effect  of  SNR.  Figure  IV-7  shows  coherence  and  phase  esti- 
mates of  a two-channel  shifted  white  signal  with  five  levels  of  independent  white  noise  added  to 
each  channel.  It  is  obvious  from  the  figure  that  the  accuracy  of  the  estimates  deteriorates  as 
uncorrelated  noise  power  increases.  The  coherence  estimate  seems  particularly  sensitive  to 
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the  presence  of  the  noise.  A comparison  with  estimates  obtained  from  a 100-sample  Bartlett 
estimator  indicates  that,  even  in  the  presence  of  noise.  Burg's  estimator  can  still  provide  more 
accurate  coherence  and  phase  estimates  than  the  conventional  estimator. 

A.  T.-U.  Ng 
T.  E.  Landers 


B.  A MULTIPLE  EVENT  AT  SEMIPALATINSK 

The  practical  application  of  seismic  discrimination  techniques  depends  to  a large  extjent  on 
both  the  natural  and  explosion  seismic  history  of  the  region  in  which  explosions  might  occur. 
Consequently,  for  a given  site  where  events  of  both  types  have  occurred  and  their  data  reduced 
to  a set  of  discrimination  parameters  (e.g.,  Mg-m^,  complexity,  3rd  moment,  etc.),  a npw  event 
can  be,  with  a great  degree  of  certainty  above  some  magnitude  threshold,  classified.  In  the  in- 
stance where  a testing  region  is  aseismic  and  new  events  are  found  to  have  the  same  parameters 
as  previous  events,  then,  through  the  somewhat  more  subtle  implication  that  earthquakes  would 
have  sufficiently  different  discrimination  parameters,  they  too  could  be  classified.  In  ar|y  case, 
when  the  depth  of  an  event  can  be  established  at  greater  than  the  maximum  feasible  testing  depth 
then  the  event  would  necessarily  be  an  earthquake,  irrespective  of  the  seismic  history  of*  the 

epicentral  region.  However,  it  is  possible  to  mimic  the  depth  phase  by  multiple  explosions  such 

1 5 

that  casual  observations  of  teleseismic  data  would  classify  the  events  as  a single  earthquake. 
Event-location  techniques,  including  master-event-location  methods,  would  be  very  difficult  to 
use  to  establish  depth  since  time  picks  on  later  phases  would  have  errors  due  to  masking  by  the 
coda  of  the  earlier  events.  Consequently,  depth  determination  for  unusually  complex  events  in 
aseismic  regions  would  depend  on  analyzing  waveform  data.  On  the  other  hand,  no  regiqn  of  the 
globe  can  be  truly  aseismic,  i.e.,  said  not  to  be  capable  of  having  an  earthquake.  In  fac^,  from 
a geologic  point  of  view,  the  earth's  crust  is  essentially  faulted  everywhere.  The  source  struc- 
ture and  stress  environment  and  path  characteristics  determine  many  of  the  teleseismic  proper- 
ties (e.g.,  Novaya  Zemlya  events  are  more  complex  than  Semipalatinsk  events,  Semipalatinsk 
events  are  much  higher  frequency  at  NORSAR  than  in  North  America,  and  NTS  events  often  re- 
lease considerable  shear  strain  energy)  such  that  the  prediction  of  earthquake  parameters  for 
current  aseismic  regions  is  very  difficult.  The  result  of  these  considerations  is  that,  ufider  the 
circumstance  of  an  event  taking  place  in  a previously  aseismic  region,  the  possibility  of*  multiple 
explosions  can  not  be  ruled  out  on  the  basis  that  its  discrimination  parameters  do  not  mfctch  a 
previously  determined  set.  In  light  of  the  above  discussion  we  present  data  on  one  such  puzzling 
event(s)  that  occurred  20  March  1976  at  Semipalatinsk,  USSR  at  50°N  77.3° E at  approximately 
4:03:39  GMT  at  0-km  G depth.16 

Figure  IV-8  shows  the  SP  vertical  seismograms  for  the  event  at  stations  with  the  b4st  focal 


sphere  distribution  that  can  be  obtained  in  the  data  base  available  to  us.  The  data  from  College 
(COL)  were  digitized  from  the  WWSSN  70-mm-film  chip,  those  from  Yellowknife  (YKC)  and 
Guaribidanur  (GBA)  are  single -channel  raw  data  recorded  on  analog  tape  and  then  run  through 
an  A-D  converter,  and  the  NORSAR  data  are  the  single -channel  04C  site  digital  recording.  The 
event  is  clearly  quite  unusual.  To  our  knowledge,  no  natural  seismic  event  near  this  magnitude 
has  occurred  within  several  hundred  kilometers  of  the  epicenter.  Yet,  what  appears  to  be  P, 
pP  at  about  7 sec  and^sP  at  about  10  sec  indicates  its  hypocenter  to  be  at  a depth  of  about  30  km, 
and  so  it  would  be  an  earthquake.  It  is  unusual,  though  perhaps  possible,  that  in  every  quadrant 
of  the  lower  focal  sphere  the  depth  phases  are  larger  than  the  initial  arrival.  That  such  behavior 
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can  be  mimicked  by  a set  of  explosions  is  illustrated  in  P'ig.  rV-9.  A fuller  explanation  of  the 
figure  is  given  by  Landers  in  Ref.  15;  however,  it  suffices  here  to  say  that  the  bottom  trace  is 
the  sum  of  the  5 upper  identical  traces  (except  for  a delay  and  positive  scaling  factor),  the  basic 
seismogram  being  a LASA  recording  of  a Semipalatinsk  explosion.  Not  only  is  there  an  appear- 
ance of  depth  phases  (at  an  earthquake-only  depth),  but  pP  is  inverted.  Figure  IV- 10  shows  the 
surface  wave  that  would  be  recorded  at  a distance  of  12°  with  the  same  scaling  and  delays  as  in 
Fig.  IV-9.  Characteristically,  signal  busting  at  the  higher  frequencies  produces  a sum  that  has 
the  overall  summed  amplitudes  at  longer  periods,  but  reduced  amplitudes  at  shorter  periods. 

The  overall  effect  is  a more  spread-out  surface-wave  train  and  a larger  Mg.  Since  the  m^  would 
still  be  obtained  from  the  first  arrival,  the  M -m^  for  the  multiple  set  would  be  pushed  toward 
earthquake-like  values.  Figure  IV- 11  shows  the  LP  data  for  the  20  March  1976  event(s)  and  a 
4 July  1976  event  with  the  same  epicenter  at  the  MAIO  SRO.  The  20  March  event  shows  the  same 
property  as  the  summed  surface  wave  in  Fig.  IV- 10,  while  the  more  compact  Rayleigh  wave  typ- 
ical of  a single  explosion  is  exhibited  by  the  4 July  event.  The  Love-to-Rayleigh-wave  ratio  in 
the  20  March  event  is  about  3:2,  whereas  for  the  4 July  event  it  is  about  1:2.  Though  more  an 
indicator  of  tectonic  conditions  at  the  epicenter,  the  variation  in  this  ratio  indicates  that  if  it  was 
a set  of  explosions,  an  extremely  large  amount  of  SH  motion  was  generated. 


TABLE  IV-1 

STATION  M VALUES  FOR  20  MARCH  1976  EVENT 
s 

Station 

A(m) 

T (sec) 

M 

s 

MAIO 

0.88 

19 

4.  17 

KON 

0.67 

18 

4.  30 

SHI 

0.67 

21 

4.  24 

TAB 

0.33 

20 

3.86 

QUE 

0.5 

22 

4.02 

KBL 

1.08 

19 

4.  17 

Mean  M = 4.  1 3 
s 

S.D.  = 0.  16 

Table  IV-1  shows  data  used  to  compute  an  M value  for  the  20  March  1976  event.  The  value 

s 

was  determined  using  the  formula  given  in  Marshall  and  Basham.  In  Fig.  IV- 12,  the  resulting 
Mg-m^  value  (point  A)  is  superimposed  on  Marshall  and  Basham's  data  events  in  the  same  epi- 
central  region.  Assuming  4 multiple  events  (for  reasons  we  will  show  when  the  cepstra  for  this 
event  are  discussed),  each  contributing  equally  to  the  surface  waves,  the  Mg-mb  typical  of  the 
individual  events  would  be  at  point  B.  As  noted  in  Fig.  rV-8,  the  secondary  arrivals  were  bigger 
than  the  primary.  Assuming  an  average  amplitude  per  arrival  to  be  1.5,  the  amplitude  of  the 
first  arrival  (the  "average"  Mg-m^  per  event)  would  fall  at  point  C,  squarely  on  the  explosion 
trend  line. 
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As  mentioned  above,  4 events  were  assumed  to  be  contributing  to  the  surface  wave;; 
dence  for  this  minimum  number  comes  from  cepstral  analysis  of  the  first  few  seconds 
SP  records.  We  note  that  this  part  of  the  seismogram  in  each  of  the  records  shown  in 
is  more  complicated  than  is  usually  observed  for  events  from  this  region.  Figure  IV- 


Of 


the  first  4 sec  of  the  Yellowknife  seismogram,  and  Fig.  IV- 14  shows  the  cepstrum  com 
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by  analyzing  the  log  spectrum  with  maximum-entropy  (ME)  spectral  analysis*""  and  by 
FFT  of  the  log  spectrum.  The  interpretation  of  this  function  is  that  2 events,  the  seco 
larger  than  the  first  and  delayed  by  0.9  sec,  are  present  and  that  depth  phases  for  both 
0.8  sec.  Such  a depth-phase  time  delay  indicates  that  the  events  originate  at  less  than 
kilometers  hypocentral  depth.  Counting  these  two,  and  assuming  that  the  later  arrivals1  are 
each  a separate  event,  makes  at  least  four  and  leads  to  the  interpretation  that  the  individual 
events  have  explosion-like  Mg-m^,  signal  characteristics,  and  depth. 

In  summary,  the  SP  data  indicate  that  more  than  2 events  occurred.  LP  data  indick 
the  events  have  more  explosion-like  Mg-m^  values  than  the  nearest  earthquakes.  The 
ity  exists  that  a combination  of  explosions  and  high-stress  drop  earthquakes  occurred  o 
series  of  extremely  high-stress  drop  very  shallow  earthquakes  occurred,  but  these  see 
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T.  E.  Landers 
M.  W.  Shields 

C.  SPLITTING  OF  THE  pP  PHASE  FROM  SUBOCEANIC  EVENTS 

The  intent  of  this  project  was  to  observe  a splitting  of  the  depth  phase  (pP)  from  suboceanic 
events.  Splitting  might  be  expected  since  energy  from  the  pP  wave  incident  on  the  oceap  floor 
should  be  partially  reflected  and  partially  transmitted.  The  part  entering  the  water  col(min 
would  be  reflected  at  the  surface,  and  would  return  to  the  ocean  floor  with  a delay  equal  to  the 
two-way  travel  time  in  the  water  column.  This  energy  would  again  be  partly  reflected  qnd  partly 
transmitted  at  the  ocean  floor,  the  transmitted  part  giving  rise  to  the  presumed  splitting  of  the 
pP  phase. 

Theoretically,  energy  would  be  trapped  in  the  water  column  by  the  strongly  reflecting  inter- 
faces at  the  ocean  surface  and  floor.  This  situation  would  produce  a periodic  repetition  of  the 
pP  phase.  The  period  would  be  the  two-way  travel  time  in  the  water  column,  and  the  amplitudes 
of  the  "multiples"  would  be  expected  to  decay  geometrically.  Removal  of  these  reverberations 
is  necessary  for  determining  the  position  of  deep  reflectors  under  the  ocean  floor.  Additionally, 
the  characteristics  of  the  reverberations  provide  information  on  the  nature  of  the  ocean  floor. 

In  searching  for  an  event  which  might  be  expected  to  exhibit  reverberations,  there  fcre  a 
number  of  criteria  which  the  event  should  meet  to  facilitate  observation  of  the  water-column 
phases.  First,  the  pP  wave  should  emerge  in  a region  under  water.  This  is  not  exactl/  a facile 
observation,  since  the  point  of  emergence  may  be  a considerable  distance  from  the  epicenter 
for  a very  deep  event.  Second,  the  event  should  have  an  uncomplicated  source  mechanism  (an 
explosion  would  be  ideal),  so  that  there  would  be  little  question  whether  observed  features  of  the 
seismogram  are  due  to  the  source  or  to  reflecting  structures  in  the  medium  of  propagation. 
Third,  it  is  desirable  that  the  epicenter  be  located  in  a region  where  the  depth  of  water  is  con- 
stant for  hundreds  of  kilometers  in  all  directions,  so  that  the  two-way  travel  time  can  be  accu- 
rately computed  even  if  the  epicenter  or  point  of  pP  emergence  is  uncertain.  This  condition 
can  be  relaxed  for  shallow  focus  events.  Fourth,  events  with  abnormally  complicated  pP  phases 
(relative  to  simpler  P phases)  may  be  of  interest.  Finally,  the  event  should  be  in  deep  water 
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to  facilitate  observation  of  the  reverberant  phases  — the  deeper  the  better,  since  the  reverbera- 
tions will  be  more  widely  separated. 

The  event  which  best  met  these  conditions  was  a shallow  event  (49  km)  in  the  Kurile  Island 
arc  (15  July  1973,  14:06:50.1,  43.1°N  x 146.52°  E,  m^  = 5.4).  A seismogram  for  this  event  from 
a single  station  at  LASA  is  shown  in  Fig.  IV-15.  Following  the  pP  phase,  this  seismogram  has 
a complicated  structure.  This  is  probably  due,  in  part,  to  arrival  of  the  sP,  PcP,  pPcP,  and 
sPcP  phases.  An  interesting  feature  of  this  seismogram  is  the  pP  arrival  itself.  Note  that  the 
P wave  has  a very  large  positive -amplitude  spike.  Presumably,  then,  the  pP  wave  should  have 
a very  large  negative-amplitude  spike.  This  pP  appears  to  have  two  such  spikes  separated  by 
about  1.9  sec.  Obviously,  some  technique  more  sophisticated  than  visual  inspection  must  be 
employed  to  determine  if  the  first  part  of  the  pP  arrival  is  actually  two  slightly  separated 
arrivals. 

Many  techniques  for  the  detection  of  multiples  could  be  applied  to  this  problem.  Those  dis- 
cussed here  by  no  means  exhaust  the  possibilities.  The  two  techniques  tried  as  part  of  this  proj- 
ect fall  into  the  class  of  cepstral  analysis  methods.  Both  involve  computation  of  the  cepstrum 

of  pertinent  parts  of  the  seismogram  and  search  for  peaks  in  the  cepstrum.  Since  details  of  the 
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theory  and  computation  ' of  the  cepstrum  are  amply  covered  elsewhere,  they  are  omitted 
here.  However,  it  is  worthwhile  to  recall  a few  properties  of  the  cepstrum  that  are  important 
to  the  present  efforts  at  interpretation. 

The  computation  of  the  cepstrum  is  not  a trivial  task,  due  to  the  fact  that  the  logarithm  of 

the  spectrum  is  a multivalued  function.  Since  computer  arctangent  routines  will  return  only  the 

principal  value  of  the  phase,  the  principal  value  must  be  ,T unwrapped"  to  produce  an  analytic 

function.  A particularly  effective  algorithm  for  performing  this  unwrapping  operation  has  been 
20 

developed  by  Tribolet. 

The  important  feature  of  the  cepstrum,  so  far  as  multiple  detection  is  concerned,  is  that 
the  cepstrum  of  a signal  with  multiples  exhibits  peaks  located  at  the  delay  times  of  the  multiples. 
If  a primary  arrival  is  larger  than  the  secondary,  then  the  cepstrum  will  be  a causal  spike  train. 
If,  instead,  the  reverse  is  true  then  the  cepstral  spike  train  would  be  anti-causal.  The  situation 
is  much  more  complicated  if  there  are  multiple  arrivals  with  different  arrival  times.  This  will 
complicate  the  determination  of  arrival  times,  although  the  peaks  at  the  individual  arrival  times 
should  be  the  largest. 

Computation  of  the  cepstrum  is  equivalent  to  spectral  analysis  of  the  log  spectrum  of  the 
signal.  One  way  of  performing  spectral  analysis,  as  mentioned  above,  would  be  to  compute  the 
Fourier  transform  of  the  log  spectrum  directly.  Harmonics  of  the  log  spectrum  should  be  well- 
resolved  in  the  cepstrum,  provided  that  the  available  bandwidth  of  the  signal  spectrum  is  large 
compared  with  the  periods  of  the  harmonics.  Of  course,  the  narrower  the  available  spectral 
bandwidth,  the  worse  the  resolution  will  be. 

In  particular,  for  a multiple  with  small  delay  time  n , only  a few  periods  of  the  fundamental 
-ju>n 

e will  be  present  in  the  log  spectrum.  The  peak  due  to  this  fundamental  may  not  be  well- 

resolved  in  the  cepstrum.  To  obtain  a better  estimate  of  the  location  (or  existence)  of  such  a 

peak,  it  may  be  necessary  to  use  one  of  the  high-resolution  spectral-estimation  techniques. 

21 

The  MEM  has  been  used  in  this  manner  by  Landers  to  determine  the  (necessarily  short)  delay 
times  of  the  pP  phases  from  underground  nuclear  explosions.  The  MEM  essentially  fits  a ratio- 
nal (denominator  only)  model  to  the  power  spectrum  of  a signal.  Since  harmonics  present  in  a 
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signal  contribute  poles  to  the  z-transform  of  the  signal,  the  MEM  tends  to  resolve  these  har- 
monics well.  It  is  known  to  provide  good  resolution  when  only  a few  cycles  of  a harmonic  are 
present  in  the  data.  The  MEM  algorithm  applied  to  the  log  spectrum  of  a signal  generates  a 
rational  approximation  to  the  squared  magnitude  of  the  cepstrum. 

Even  the  best-quality  digital  seismic  data  are  rather  poor  for  the  detection  of  closely  spaced 
arrivals,  because  the  data  are  very  narrowband.  For  LASA  SP  data,  the  available  bandwidth  is 
approximately  2 Hz,  even  though  sampling  is  at  20  Hz.  A multiple  arriving  with  a delay  of  2 sec 
would  produce  a harmonic  in  the  spectrum  with  only  4 cycles  in  the  available  spectrum.  This 
is  the  reason  that  high-resolution  spectral  estimation  must  be  employed  to  estimate  the  delay. 

The  processing  of  the  Kurile  Island  event  is  summarized  in  Figs.  IV- 16  through  IV-2  0. 
Cepstra  computed  directly  using  Fourier  transforms  are  shown  in  Figs.  IV-16  and  IV-17j  Some 
filtering  of  the  log  spectrum  was  performed  to  remove  the  low-time  portion  of  the  cepsttfum. 

It  is  probably  dominated  by  the  source  mechanism  of  the  earthquake.  Some  of  the  high-tlme 
part  of  the  cepstrum  was  removed  also.  Less  was  removed  in  the  second  example  (Fig.  (IV- 17) 
than  in  the  first.  In  the  computation  of  the  cepstra,  only  the  log  magnitude  of  the  spectrum  was 
used.  This  accounts  for  the  fact  that  the  cepstra  are  even.  The  cepstra  are  also  rather!  noisy, 
which  may  be  due  in  part  to  the  fact  that  the  data  were  not  modulated  to  close  the  gap  around 
zero  frequency. 

Both  cepstra  are  dominated  by  large  negative  spikes  at  10.7  sec,  which  are  undoubtedly  due 
to  the  pP  arrival.  There  are  also  negative  spikes  around  12.4  and  12.6  sec;  these  are  consistent 
with  the  first  water  reverberation  phase.  In  Fig.  IV-17,  the  spike  at  34  sec  may  be  due  |o  PcP. 

Figures  IV-18  through  IV-20  show  the  steps  in  the  computation  of  MEM  cepstra.  The  log 
magnitude  and  phase  of  the  spectrum  are  highpass  filtered  to  remove  source-mechanism  effects. 
Then  the  log  spectrum  is  windowed  to  retain  the  portion  with  significant  energy  and  discard  the 
remainder.  Any  linear  trend  is  removed  from  the  windowed  portion  before  an  MEM  estimate 
of  the  cepstrum  is  computed.  Rational  approximations  of  orders  60  and  99  are  presented  in 
Figs.  IV-18  and  IV-19,  respectively.  In  Fig.  IV-20,  the  seismogram  is  windowed  to  isolate  the 
pP  arrival  and  its  coda,  and  the  cepstral  estimate  is  computed  for  pP  and  its  coda. 

From  the  Herrin  tables,  this  event  should  have  a PcP  arrival  at  23  sec.  sP  should  arrive 
around  5 sec  after  pP,  pPcP  should  be  34  sec  after  P,  and  sPcP  should  be  38  sec  after  P.  In 
Fig.  IV-18,  there  is  an  arrival  at  12.5  sec;  this  may  be  an  amalgam  of  pP,  sP,  and  the  water 
phase.  The  spike  at  2 3 sec  is  probably  PcP.  The  spike  at  5 sec  corresponds  to  pP-sP. 

Figure  IV-19  has  small  peaks  at  34  and  38  sec;  these  could  be  due  to  pPcP  and  sPcP.  The 
spike  at  12  sec  is  probably  pP.  The  large  spike  at  5 sec  is  again  consistent  with  pP-sP. 

Figure  IV-20  has  a pulse  at  —2  sec;  this  could  be  the  water  phase.  The  peak  at  4 s£c  is 
probably  sP-pP,  in  which  case  the  arrival  at  8 sec  would  be  the  second  harmonic.  Its  amplitude 
is  correct  for  this.  The  peaks  at  larger  times  are  in  the  PcP  regions. 

The  case  for  the  water  reverberation  phase  seems  to  be  inconclusive,  although  somp  of  the 
other  arrivals  may  contribute  peaks  to  the  cepstra.  If  it  is  present,  then  it  most  probably  is 
delayed  by  about  2 sec  with  respect  to  pP  as  is  suggested  by  the  raw  seismogram.  The  epicenter 
for  this  event  is  located  under  1000  to  2000  m of  water,  and  a two-way  travel  time  of  2 dec  trans- 
lates into  a depth  of  1500  m.  That  delay  is,  at  least,  consistent. 

Since  the  P and  pP  \yaves  are  well-separated  in  this  event,  it  would  be  worthwhile  to  try  to 
design  a Wiener  spiking  filter  using  the  P wave  and  then  to  apply  it  to  the  pP  wave  and  coda. 

One  further  experiment  that  might  be  tried  would  be  to  subtract  the  cepstrum  of  the  P wave  from 
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the  cepstrum  of  the  pP  wave  and  coda.  In  theory,  that  should  eliminate  the  source-function  ef- 
fects and  leave  just  the  reflector  information  in  the  pP  cepstrum.  Certainly,  more  examples 
and  a wider  variety  of  techniques  should  be  used  to  establish  the  existence  or  nonexistence  of 

water  phases.  D.  B.  Harrist 

T.  E.  Landers 


D.  POLARIZATION  FILTERING  OF  LONG-PERIOD  SRO  DATA 

As  part  of  an  evaluation  of  SRO  data,  we  are  investigating  a wide  dynamic  range  of  body  and 
surface  waves  recorded  at  Mashad  (MAIO)  from  an  aftershock  sequence  in  the  Kuriles.  These 
d^ta  are  described  with  examples  in  Sec.  Ill  of  this  SATS. 

One  promising  type  of  analysis  is  polarization  filtering  which  has  previously  been  applied 
22  23 

to  SP  data.  * In  the  SRO  network  it  is  important  to  take  advantage  of  the  3 -component,  LP 
data  at  each  station  to  enhance  body-  and  surface-wave  arrivals.  Polarization  filtering  is  par- 
ticularly useful  for  this  because  it  allows  one  to  calculate  azimuth  and  ellipticity  of  detected 
surface-wave  trains,  and  the  azimuth  and  dip  angle  of  enhanced  body  phases.  Thus,  polariza- 
tion filtering  of  3 -component  data  at  a single  station  can  reveal  some  of  the  seismic  parameters 
usually  obtained  by  array  analysis. 

Polarization  filtering  in  the  time  domain  is  based  on  the  geometry  of  the  least-squares 
particle-motion  ellipsoid  shown  in  Fig.  IV-21,  calculated  from  a window  of  3-component  data. 
The  direction  and  squared  length  of  each  axis  are  given  by  an  eigenvector  and  the  corresponding 
eigenvalue  of  the  3x3  covariance  matrix  R = {r^}  where 


r. . 


n 

I xi(t)  Xj(l) 
t = l 


(IV-1) 


In  this  equation,  n is  the  number  of  points  in  the  window;  and  x^  (t) , x^ft),  and  x^(t)  are,  respec- 
tively, the  seismogram  components  Z,  R,  and  T shown  in  Fig.  IV-21.  As  the  window  slides  in 
time  through  the  data,  the  direction  and  magnitude  of  the  largest  axis  i called  the  principal 
component  of  the  data,  constantly  change.  If  k ^ at  a given  time,  then  the  data  are  strongly 

polarized  along  the  H ^ direction.  Examples  of  linearly  polarized  phases  are  P,  S,  and  Love 
waves. 
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A standard  method  * of  enhancing  strongly  polarized  data  is  to  multiply  the  data  compo- 
nents at  the  center  of  the  window  by  a gain  factor 

g=(i-i1/i2)  . o«g<i  (rv-2) 


and  then  project  the  principal-component  vector  onto  these  scaled  data  components.  This  method 
clearly  enhances  strong  body  waves,  but  in  many  cases  it  also  enhances  noise  and  surface  waves 
which  have  ellipticities  significantly  larger  or  smaller  than  1,  and  which  are  uncorrelated  on  2 
or  more  components. 

Figure  IV-22  shows  an  example  of  LP,  3 -component  data  for  an  earthquake  in  the  Kurile  Is- 
land sequence  recorded  at  MAIO.  The  earthquake  had  5.4,  origin  time  8:07:10  on  22  Janu- 
ary 1976,  with  a distance  A = 66°  from  MAIO.  The  first  three  traces  are  the  unfiltered  Z,  R, 
and  T components,  and  the  next  three  traces  are  bandpass-filtered  data  components  (11.11  to 
33.33  sec)  using  a phase-free  operator. 

t Department  of  Electrical  Engineering,  M.I.T. 


86 


gh  waves, 
P wave 


The  3 -component  data  were  filtered  to  pass  linearly  polarized  phases  using  a sliding  win- 
dow 30  sec  wide.  The  results  are  shown  in  traces  8,  9,  and  10.  Trace  7 shows  the  filler  gain 
g given  by  Eq.  (IV-2)  as  a function  of  time.  As  the  window  moves  by  the  P,  S,  and  SS  pjhases 
the  gain  approaches  1.0  and  the  principal  component  is  passed  and  projected  onto  the  Z„  R,  and 
T axes.  Later  in  time  one  can  see  the  enhancement  of  Love  waves  on  the  T component  and  a 
reduction  of  the  Rayleigh-wave  amplitude  by  the  polarization  filtering.  The  Rayleigh  waves  are 

attenuated  because  i . and  f 9 are  nearly  equal  so  that  gain  g is  close  to  zero.  However!  the  SS 

1 L 24 

wave  train  which  contains  body  phase  and  leaky-mode  surface  waves  is  passed  because  the 

R component  is  much  larger  than  the  Z component,  causing  g to  be  clolse  to  1. 

An  interesting  problem  is  that  of  enhancing  elliptically  polarized  surface  waves  while  re- 
jecting body-wave  motion.  This  is  not  possible  using  only  the  particle^-motion  ellipsoid,.  One 
new  modification  of  the  time -domain  polarization  scheme  appears  to  be  useful  for  this  problem. 
This  is  to  reduce  the  components  to  Z and  R only,  Hilbert  transform  the  R components,  and 
use  the  correlation  coefficient  between  components  as  a gain  factor.  For  P and  Raylei 
we  expect  particle-motion  diagrams  similar  to  those  displayed  in  Fig.  IV-23(a-b).  The 
is  linearly  polarized,  but  the  Rayleigh  wave  is  elliptically  polarized  with  an  ellipticity  vVhich 
strongly  depends  on  its  propagation  mode  and  the  crustal  structure  of  the  travel  path.  There- 
fore, no  simple  relation  between  t ^ and  i ^ exists  which  is  useful  for  enhancing  surface  waves. 

However,  for  flat-layered  isotropic  crustal  models,  the  R and  Z conponents  are  90°  o*ut  of 
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phase  for  all  frequencies.  If  the  R component  is  Hilbert  transformed,  it  should  be  in  phase 
with  the  Z component  over  the  whole  dispersed  wave,  and  the  resulting  particle  motion}  will  then 
be  linearly  polarized  with  i ^ » i^.  Conversely,  body  phases  which  were  in  phase  shoi  ld  be  90° 
out  of  phase  and  appear  to  be  elliptically  polarized. 

This  is  illustrated  in  Fig.  IV-24  using  the  Z and  R components  of  Fig.  IV-22.  Trance  4 is 
the  R component  of  trace  2 after  Hilbert  transformation.  Traces  3 and  4 show  P and 
which  are  90°  out  of  phase,  surface-wave  components  in  SS  which  are  180°  out  of  phase, 

Rayleigh  phases  which  are  in  phase. 

Traces  6 and  7 in  Fig.  IV-24  show  the  results  of  polarization  filtering  traces  3 and  4 using 
the  ellipticity  gain  g in  Eq.  (IV-2)  to  pass  surface  waves.  The  R component  in  trace  7 has  been 
inverse  Hilbert  transformed  to  restore  the  original  phase  in  trace  2.  Traces  8 and  9 show  the 
results  of  subtracting  the  passed  surface  waves  in  traces  6 and  7 from  the  original  data.  Although 
the  Rayleigh  waves  are  removed,  the  results  for  body  phases  are  not  convincing.  In  particular, 
the  P phase  in  the  Z component  of  trace  6 is  not  attenuated  much,  and  on  traces  8 and 
P wave  has  only  about  one-half  the  amplitude  of  the  unfiltered  data. 

A further  refinement  in  the  polarization  scheme  is  to  replace  the  gain  factor  in  Eq.  (IV-2)  by 
the  correlation  coefficient 


S phases 
and 


g = 


12 
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9 the 


(IV-3) 


where  the  covariance  r.j  is  given  by  Eq.  (IV- 1).  This  has  the  effect  of  attenuating  components 
which  may  be  strong  but  uncorrelated,  such  as  certain  types  of  noise  or  components  which  are 
90°  out  of  phase.  Using  this  gain-factor  polarization  filtering  of  traces  3 and  4 to  pass  surface 
waves  yields  traces  11  and  12,  and  traces  13  and  14  show  the  original  data  with  surface  waves 
subtracted.  Traces  11  and  12  show  almost  complete  attenuation  of  the  P phase,  because  of  the 
90°  phase  shift  between  R and  Z components  in  traces  3 and  4.  Also,  the  long  wave  train  asso- 
ciated with  the  SS  arrival  clearly  separates  into  surface-wave  and  body-phase  arrivals  on  the 
two  pairs  of  polarized  components. 


This  combination  of  using  the  principal  component  of  the  particular  motion  ellipse  and  a 
gain  factor  given  by  the  correlation  coefficient  between  components  appears  to  produce  better 
separation  of  body-  and  surface-wave  phases  than  is  possible  using  only  the  least-squares  el- 
lipse parameters. 

This  technique  will  be  applied  to  the  problems  of  separating  interfering  phases  of  the  Kurile 
sequence  and  lowering  the  detection  threshold  of  body  and  surface  waves. 

C.  W.  Frasier 
R.  G.  North 


E.  AN  EXACT  SOLUTION  TO  THE  PROBLEM  OF  EXCITATION 
OF  NORMAL  MODES  BY  A PROPAGATING  FAULT 

Consideration  of  both  spatial  and  temporal  dimensions  of  a seismic  source  as  well  as  the 
direction  of  rupture  propagation  becomes  a significant  factor  in  excitation  of  normal  modes  and 
surface  waves  when  the  wavelengths  or  frequencies  of  the  waves  become  greater  than  about  1/10 
of  either  the  source  length  or  the  source  duration. 

Ben-Menahem2^  presented  a solution  to  the  problem  of  excitation  of  surface  waves  by  prop- 
agating sources,  valid  when  the  distance  to  the  receiver  is  large  in  comparison  with  both  the 
wavelength  and  the  source  dimension.  This  assumption  clearly  will  not  be  satisfied  in  many 
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cases  of  interest  in  seismology.  Kanamori  and  Anderson  circumvented  this  difficulty  by  di- 
viding the  800-km-long  fault  of  the  i960  Chilean  earthquake  into  16  segments,  each  represented 
by  a point  source  activated  with  a proper  time  delay.  Although  such  a procedure  leads  to  satis- 
factory results  within  the  range  of  wavelengths  considered  in  their  study,  it  is  very  time  con- 
suming if  applied  to  a great  many  modes  and  stations.  Clearly,  an  analytical  solution  is  desir- 
able. Even  if  such  a solution  should  prove  to  be  not  economical  in  practice,  it  would  allow 
testing  the  precision  and  limits  of  applicability  of  various  approximations. 

Here,  we  shall  present  an  exact  solution  to  the  problem  of  excitation  of  normal  modes  by  a 

fault  propagating  with  a constant  velocity  along  a fraction  of  a great  circle. 
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Gilbert  and  Dziewonski,  following  the  theoretical  development  of  Gilbert,  gave  the  equa- 
tions for  excitation  of  normal  modes  by  an  arbitrary  seismic  source.  The  displacement  vector 
t h 

u for  a k normal  mode  observed  at  location  r is,  in  the  time  domain: 

A' 

i 

uk(r,  t)  = ? s ™(r)  *™(t)  * Ck(t)  (IV-4) 

m =-f 

or  in  the  frequency  domain: 

i 

uR(r,  w)  = Ck(u>)  £ s™(r)  * ™(w)  (IV- 5) 

m = -i 

where  i is  the  angular  order  number  of  the  mode.  The  advantage  of  this  representation  is  that 
the  receiver  term  ^(r)  and  the  source  term  are  clearly  separated;  the  source  term  remains  the 
same  if,  for  example,  the  strain  tensor  at  ^r  instead  of  the  displacement  vector  were  to  be 
evaluated. 

The  source  function  is: 

dv  M(_Tq,  t)  : Ekm(r0)  (IV-6) 

"V0 
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where  E^  (rQ)  = l/2(Vs^  + s^  V);  the  bar  denotes  a complex  conjugate  and  M(rQ,  t)  is  t(ie  mo- 

ment density  rate  tensor. 

Gilbert  and  Dziewonski  concentrated  their  attention  on  point  sources  (in  space): 

M(r0.t)  = M(t)  • 6(r-  rQ) 

and  Eq.  (IV-6)  simplified  to 

* “(t)  = M(t)  : E™(r0)  . 

For  large  earthquakes  with  fault  lengths  of  the  order  of  several  hundred  kilometers,  it  ife  im- 
portant to  account  for  the  finite  dimensions  of  the  source  as  well  as  the  finite  time  and  direction 
of  rupture  propagation.  The  fault  length  is  clearly  the  dominant  dimension  for  such  an  darth- 
quake,  and  it  is  reasonable  to  approximate  the  volume  integral  Eq.  (IV-6)  by  a line  integral  eval- 
uated along  the  fault  of  a length  L: 


= £ j’0L  » *)  : Ekm(4)  d| 


(rv-7) 


We  shall  further  assume  that  the  equivalent  source  radius  (depth)  remains  constant  (rQ),  and  that 
the  source  line  can  be  described  by  the  arc  of  a great  circle  with  the  end  points  having  CDlati- 
tudes  and  longitudes  (0q,  anc*  resPectively*  For  a spherically  symmetric,  nonrotat- 

ing earth  the  choice  of  a coordinate  system  is  arbitrary,  and  we  may  choose  the  one  that  will 

30 

make  evaluation  of  Eq.  ( IV —7 ) the  most  convenient.  Following  Backus,  we  rotate  the  coordinate 
system  such  that  0Q  and  AQ  correspond  to  a colatitude  of  90°  and  longitude  of  0°  in  the  ndw  co- 
ordinate system,  which  we  shall  call  "equatorial,"  and  (0^,  A^)  are  (90°,  <i>),  where  4>  = L/r^. 

Note  that  the  moment  rate  tensor  must  also  be  represented  in  this  new  coordinate  systerp.  Thus, 
the  line  source  represents  now  a fraction  of  the  equator.  Therefore, 


Ot) 


I Mto.t)  : E™(r0>£,<p)  dtp 


For  a uniform  slip  motion  with  a constant  propagation  velocity  v: 

<Prnx  / (Prc 


also: 


Thus, 


/ ^ 0\  / ^ 0 \ 
M(<p,t)=  M(0,t—  — M0(t-  -£)  ; 


7r  N 7m,  7r  % , -incup 
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1 r n/i  u <pro\  ~m/  ,r\ 

4>  J0  ~0  V V / : ~k  (r0’  2 ) 


e-im^  d«p 


(IV-8) 


and  in  the  frequency  domain: 


-=  J M0(w)  : 6™(rQ,  |)  exp{-i(p  [(ar0/v)  + m]}  dtp 


M / x . 7 m,  7c_ v sln  *m  "lxm 

= ’ £k  (r0>  2 v e 
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(IV— 9) 
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whe  re 


*m=  f*(-F  +m)  * <IV-10) 

Function  of  Eq.  (IV -9)  may  now  be  substituted  into  Eq.  (IV-5)  and  yield  the  desired  solu- 

tion to  the  problem. 

The  "directivity"  term  dependent  on  x is  similar  in  its  form  to  that  given  by  Ben- 
2 6 

Menahem  for  the  excitation  of  surface  waves  by  propagating  sources.  Following  some  ele- 
mentary transformations,  xm  can  be  written  as: 


m 


u>L  / r0 

= 2C^[  — 
r0  ' 


m 


i + 1/2 


where  = u>r0/(i  + l/2)  is  the  phase  velocity  at  a radius  rQ.  For  shallow  sources,  Cr  » C, 

the  phase  velocity  at  the  surface  of  the  earth. 

Ben-Menahem’ s directivity  term  is 

~ L/C  v 

XR  “ 2C  ( v COSI/) 


where  v is  the  angle  measured  at  the  origin  of  the  rupture  between  the  direction  of  rupture  prop- 
agation and  direction  of  the  receiver. 

The  equivalence  between  cos  v and—  m/(i  + l/2)  can  be  established  on  the  basis  of  the  as- 
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ymptotic  properties  of  spherical  harmonics.  Dahlen  has  shown  that  for  a high  angular  order 

number  i,  a spherical  harmonic  ¥^(0,  (p)  can  be  associated  with  a directed  great-circle  path 

with  a pole  at  (0,  A),  where  cos0  = m/(i  + l/2).  There  is  a 180°  difference  between  the  angle  0 

defined  by  Dahlen  and  that  in  our  coordinate  system. 

For  evaluation  of  amplitudes  of  singlets  in  the  geographical  coordinate  system,  we  must 
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perform  the  back -transformation  to  the  original  source  coordinates.  Backus  gives  recurrence 
formulas  that  provide  a very  efficient  means  of  computation  of  the  elements  of  the  rotation  ma- 
trix. If  one  desires  to  obtain  amplitudes  of  the  degenerate  multiplets,  the  equatorial  coordinate 
system  is  the  most  efficient  for  such  computations. 

A.  M.  Dziewonski 

B.  A.  Romano  wicz 

F.  DISSIPATION  OF  SHEAR  AND  COMPRESSIONAL  ENERGY  IN  THE  MANTLE 
FROM  NORMAL-MODE  DATA 

Observations  of  the  attenuation  of  modes  of  free  oscillation  provide  the  best  measure  of  the 
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average  anelastic  properties  of  the  earth.  In  an  earlier  SATS,  we  described  a technique  for 
stacking  spectra  from  many  stations  which  enabled  us  to  measure  Q for  over  200  modes.  We 
also  presented  a result  from  inversion  of  these  data  — a two-layer  model  of  the  mantle  with 
Q = 120  in  the  upper  mantle  and  Q = 240  below  670  km  depth.  We  remarked  that  this  model, 

[j.  p. 

which  must  represent  a lower  bound  for  Q , does  indeed  indicate  a much  lower  value  of  Q in 

M*  M- 

the  lower  mantle  than  any  others  in  the  literature.  For  example,  models  "MM8"  (Ref.  33)  and 

"LMS"  (Ref.  34)  have  Q = 7 50  in  the  lower  mantle. 

M’ 

Here,  we  use  new  data  to  study  the  variation  of  Q in  the  mantle.  By  inversion  experiments 
using  Q observations  for  radial  modes,  LP  spheroidal  modes,  and  high-Q  overtones,  we 
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conclude  that  the  value  of  Q in  the  lower  mantle  is 


400.  We  also  find  that  a valut 
Qgcs  = is  consistent  with  Q models  which  satisfy  free  oscillation  data  [published  meat; 
ments  of  ^5gcg  range  from  600  (Ref.  35)  to  162  (Ref.  36)].  In  what  follows,  we  present  only 
major  conclusions;  the  details  will  be  published  elsewhere. 

Q measurements  for  35  modes  were  obtained  using  individual  LaCoste  gravimeter  r0 
of  three  major  earthquakes  (1964  Alaska,  1970  Colombia,  1975  Solomon  Islands).  We  us 


cordings 
<?d  the 
oh  dif- 


simple  method  of  observing  the  time  rate  of  decay  of  a spectral  peak  corresponding  to  ea< 
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ferent  mode. 

Figures  IV- 2 5 and  IV- 2 6 show  the  application  of  this  method  to  the  m6des  qSq  and  QSk  using 

recordings  of  the  1964  Alaskan  earthquake.  The  data  consist  of  two  time  series,  recorded  at 

the  rate  of  one  sample  per  minute  on  two  LaCoste  gravimeters  located  at  UCLA  (designated  as 
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#4  and  #7).  The  series  are  over  18  days  long  and  thus  provide  an  excellent  opportunity!  for 
observing  high-Q,  LP  modes. 

The  mode  qSq  is  remarkable  for  its  very  high  Q.  Previous  estimates  of  this  value  have 
ranged  from  900  to  >25,000.  We  believe  that  the  proper  value  of  Q is  actually  less  than  7000. 
As  shown  in  Fig.  FV-25,  for  gravimeter  #4  we  obtain  the  value  5096  (467  5 to  5600),  and  fqr  gra- 
vimeter #7  we  obtain  5663  (5203  to  6212). 


+ 1 - 1 

The  value  of  Q for  qS2  is  about  550.  The  singlets  qS2  and  qS2  can  be  resolved  ai 


Q measured  separately.  The  m = 0 singlet  was  not  observable  because  of  the  geographic 
tion  of  this  station.  In  Fig.  IV-26,  the  open  circles  are  values  of  Pq,  the  maximum  power 
peak  QS2  (T  = 3287.76  sec).  Overlapping  time  segments  of  6-day  lengths  were  used.  Th£ 
line  shows  the  least-squares  fit  for  all  9 points  shown,  and  gives  Q = 581  (509  to  678  witlji: 
S.D.  of  Q"1).  If  only  the  first  7 points  are  used,  the  result  is  Q = 499  (447  to  565),  and 
4 points  alone  give  Q = 531  (520  to  542). 


__  + 1 

The  open  triangles  in  Fig.  IV-26  are  P values  for  QS2  . P is  the  average  power,  fo 


integrating  from  period  T^  to  T2  or  ±3  discrete  frequency  points  from  the  central  peak  v; 
The  line  fitted  to  the  first  7 of  these  points  gives  Q = 586  (542  to  638). 

The  solid  triangles  in  Fig.  IV-26  are  P values  from  integrating  across  the  entire  m 
(from  T^  = 3304.34  to  T2  = 3130.70,  which  corresponds  to  33  discrete  frequency  points) 
uncertainty  in  this  fit  represents  about  10  percent  of  the  measured  value  of  Q of  5 56. 

Having  thus  obtained  Q measurements  for  35  modes  using  individual  recordings  of  3 
quakes,  we  performed  inversion  experiments  to  find  Q models  which  satisfy  these  data, 
the  purposes  of  inversion,  we  added  to  this  data  set  5 modes  — 3 toroidal  and  2 spheroidal 
whose  Q values  had  been  measured  by  stacking.  We  then  had  40  modes:  21  were  fr<>: 
to  qS^5;  7 radial  modes  qSq  to  ^Sq;  8 high-Q  spheroidal  overtones;  3 toroidal  modes; 

1 Stonely  branch  mode, 

We  employed  two  methods  of  inversion  — inversion  in  the  data  space  and  in  the  pararfi 
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space.  The  method  of  Gilbert,  called  "ranking  and  winnowing,”  takes  into  account  the 
the  observational  errors  in  finding  Lagrange  multipliers  for  inversion  in  the  data  space 
these  inversion  results  is  shown  in  Fig.  IV-27.  Model  MM8  satisfies  the  data  set  of  40  m< 
with  a relative  rms  residual  of  0.416,  and  it  predicts  Q values  higher  than  all  those  obs 
The  model  obtained  by  ranking  and  winnowing  using  these  data  and  MM8  as  a starting 
has  a relative  rms  of  0.196.  Figure  IV-27  shows  that  this  model  has  < 400  in  the  low£ 
tie.  Our  inversions  in  the  data  space,  which  are  free  of  prejudice  regarding  an  initial  c 
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parameterization  of  the  model,  failed  to  reveal  any  pronounced  resolvable  features  in  the  lower 
mantle. 

The  parameter  space  inversion  was  performed  for  2 or  3 shell  models  of  Q in  the  mantle. 

For  example,  the  simple  model  Q = 353  below  670  km,  Q = 110  in  the  upper  mantle,  and  Q,v  = °o 

M-  H'  k 

everywhere  satisfies  our  data  set  with  a relative  rms  of  0.203.  However,  this  model,  like  all 

models  we  obtain  when  is  assumed  to  be  infinity,  predicts  Q values  too  high  for  all  the  ra- 
dial modes.  Thus,  we  allowed  QK  to  vary  in  subsequent  inversions.  Because  of  the  small  size 
of  this  data  set,  we  could  not  allow  both  Q and  QK  to  vary  in  more  than  one  layer  without  obtain- 
ing  a negative  value  for  QK  in  the  lower  mantle.  However,  by  assuming  that  = <*>  in  the  lower 
mantle,  but  allowing  to  assume  a finite  value  in  the  upper  mantle,  we  obtain  a much  better 
fit  to  the  data.  The  model  with  = 415  in  the  lower  mantle  and  112  in  the  upper  mantle,  with 
Qk  = 583  in  the  upper  mantle,  has  an  rms  of  0.185.  This  model  satisfies  the  observations  of  ra- 
dial modes  better  than  models  with  QK  = <*>,  and  provides  indication  that  QK  may  be  finite  in  the 
upper  mantle.  Further  inversion  studies  and  consideration  of  the  resolution  of  the  data  will  be 
needed  to  confirm  this  result.  ^ ^ Sailort 

A.  M.  Dziewonski 


G.  VELOCITY-DISPERSION  KERNELS  FOR  FREQUENCY-DEPENDENT  Q 

Recently,  numerous  authors  have  drawn  attention  to  the  problem  of  velocity  dispersion  due 
to  anelasticity.  If  the  attenuating  process  has  properties  of  a minimum  phase  filter,  then  the 
real  and  imaginary  parts  of  the  complex  wavenumber  k(u>)  are  interrelated  through  the  Kramers- 
Kronig  relationships: 

Re[nM-n(0)]=  hL  P f 4?  bnnM  (IV-lla) 

* '> 0 <7(<X  - a,  ) 

Imn(ui)  = - — P f dg  Re  n(g>  (IV-llb) 

* J0  <j  — io 


where  P denotes  the  Cauchy  principal  value;  n(u>)  = k(u>)/(u>/c),  c being  the  nondispersive  phase 
velocity;  Imn(w)  corresponds  to  l/2Q(u>)  (see  Ref.  40). 

In  inversion  studies  where  both  elastic  and  anelastic  properties  of  the  medium  are  investi- 
gated (see  Sec.  H below),  it  is  important  to  establish  the  functional  form  of  Imn(u>),  in  order  to 

evaluate  the  integral  Eq.  (IV-lla).  In  most  cases  it  is  assumed  that  Q is  frequency  independent, 
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at  least  within  the  seismic  frequency  band,  and  Eq.  (IV-lla)  then  has  a simple  solution.  * 

There  is  a growing  body  of  evidence  that  this  assumption  may  represent  a fairly  good  approx- 
imation of  reality.32,3^*42  Yet,  as  the  quality  of  the  Q data  improves,  it  may  be  necessary  to 
consider  Q that  varies  both  with  radius  and  frequency.  We  propose  here  a model  that  could  be 
useful  in  the  parameter  space  inversion  studies. 

_ A 

Let  us  define  q^  = Q.  at  a number  of  discrete  frequencies  uk,  and  let  q decrease  to  zero, 
as  a linear  function  of  the  logarithm  of  frequency,  at  and  Obviously,  q^  = q^  = 0.  In 

this  way,  q(w)  is  defined  everywhere  as  a continuous  function  of  frequency.  For  ok  4 w ^i+l* 
we  have: 


q(w) 


[q.  Jn(w/w.+1)  - qi+l  in  (w /w. )] 
in(w./u>.+1) 


t Department  of  Geological  Sciences,  Harvard  University,  Cambridge,  MA  02138. 
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The  preceding  equation  could  represent  a basis  for  a linear  inverse  problem  of  retrieving  q.,(r) 
from  attenuation  data.  The  overall  change  in  the  refraction  index  is  also  linear  with  respect 
to  q.: 

N-l 

Re  [n(o>)  — n(0)]  = | £ q.  • IjM  . 


I IV- 12) 


i=2 


Thus,  He  [n(a>)  — n(0)]  is  a linear  function  of  q^  multiplied  by  frequency-dependent  kernel^  I.(u>). 
These  kernels  are  obtained  by  evaluation  of  the  integrals: 


- 1 r 

in(«i+i/«i) 


w.  (ina  — into.^)  da 


/ 2 2. 

a(a  - u)  ) 


wi+l  (in<T"  inwi+i)  d<T 

71  2. 

a(a  — (jj  ) 


rv-i3) 


The  result  of  integration  is: 


where 


and 


in(wi/wi_1)  + *n<wi+1/wi) 


Si"Si-l 


S.  - S.i4 
1 1+1 


V 


-F(x) 


-ir2/24 


x = lj/w  . for  u)  < cj  . 

J 3 


for  u)  = 


J o ^ 

F(x)  + y (inx)  — JT  /l2  ; x = c o./u>  for  a>  > a;. 


(IV- 14) 


F(x)  = E <2^ 


1 v2  2n 
x 


n=  1 


(kv- 


15a) 


The  infinite  sum  above  is  rapidly  converging  for  x<  l/2;  for  l/2  < x < 1,  it  is  more  economical 
to  use  the  expression  (due  to  F.  Gilbert): 


2xn 

(1  - x ) 


n=  1 


IV -15b) 


In  Fig.  IV-28,  a "triangular”  q(co)  function  is  compared  with  an  absorption-dispersion  pair 
that  leads  to  the  same  overall  change  in  velocity:  v(°°)  — v(0).  Functions  that  correspond  to  a 
change  in  velocity  are  very  close  to  each  other  over  the  entire  frequency  range.  Even  the  atten- 
uation functions  show  great  similarity.  It  appears  that,  for  practical  purposes,  the  triangular 
function  q(u>)  can  be  used  as  a simulator  of  absorption-dispersion  pairs.  Their  advantage  is  that 
they  are  much  more  convenient  in  designing  a desired  Q(u>)  function. 
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In  Fig.  IV-29,  we  compare  four  different  functions  Q(u>)  and  the  relative  changes  in  velocity 

41 

associated  with  them;  following  Liu  et  aL,  we  determine 
*v/v=  . 

It  is  important  to  notice  that  the  velocity  dispersion  for  Q that  varies  with  frequency  can 
be  very  well  matched,  in  a period  range  from  1 to  10,000  sec,  by  velocity  dispersion  for  a con- 
stant Q in  the  seismic  frequency  band.  Thus,  there  is  an  equivalent  constant  Q that  gives  ap- 
proximately the  same  velocity  dispersion  as  that  computed  for  a Q(a>)  function  that  varies  by  a 

43 

factor  3 to  5 in  the  seismic  frequency  range.  Kanamori  and  Anderson  indicate  that  this  is  an 
upper  limit  on  the  range  of  variation  of  Q in  the  seismic  frequency  band. 

In  both  cases  shown,  this  equivalent  Q is  only  20 -percent  higher  than  the  average  Q in  the 
normal  mode  period  range  (100  to  3000  sec).  We  conclude  that  assumption  of  a constant  Q model 
in  inversion  studies  in  which  the  attenuation  data  are  represented  by  measurements  of  Q of  nor- 
mal modes  should  not  lead  to  grossly  erroneous  results. 

A.  M.  Dziewonski 


H.  FINITE-STRAIN  EARTH  MODEL  WITH  CONSIDERATION 
OF  VELOCITY  DISPERSION  DUE  TO  ANELASTICITY 

41  44  45 

It  has  been  pointed  out  in  several  recent  studies  * * that  the  effect  of  velocity  disper- 

sion due  to  attenuation  within  the  range  of  periods  from  1 to  3000  sec  is  several  times  greater 
than  the  errors  in  measurements  of  travel  times  of  body  waves  or  periods  of  free  oscillations. 

In  the  first  inversion  study  in  which  the  effect  of  velocity  dispersion  was  taken  into  accouht, 

46 

Hart  et  aL  used  only  toroidal  modes.  Also,  they  have  adopted  a Q model  MM8  of  Anderson 

33  ^ 

et  aL  in  which  Q in  the  lower  mantle  is  approximately  two  times  greater  from  recent  esti- 
mates based  on  new,  more  accurate  measurements  of  Q of  normal  modes. 

This  SATS  contains  preliminary  results  of  an  attempt  at  a simultaneous  inversion  for  both 
elastic  and  anelastic  parameters  of  the  earth1  s structure.  The  data  set  includes  40  high-quality 

measurements  of  Q of  normal  modes  (see  Sec.  F above),  the  free  oscillation  data  of  Gilbert 

28  47  48  49  50  51  49 

and  Dziewonski,  travel  times  for  P ’ and  S * (surface  and  deep  focus),  PKP,  SKS, 

PcP-P,52  and  PKiKP-PcP.53 

-1 

Assuming  that  q(r),  where  q(r)  = Q (r),  is  independent  of  frequency  (see  Sec.  G above)  and 
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that  it  is  our  intention  to  derive  a model  valid  at  a 1-sec  period,  a perturbation  in  a period  of 
free  oscillation  or  a travel  time  of  a body  wave  can  be  expressed  as: 


dr  (6pR  + 6V  P + 6V  S-  ^ 6q  M - 6q  K) 

P S 7T  nfJL  7T 


(IV- 16) 


where  R,  P,  S,  M,  and  K are  differential  kernels  for  perturbations  in  density,  compressional 
and  shear  velocity,  q and  q respectively;  r designates  the  period  of  the  wave  under  considera- 

jj.  * 

tion  (t  = T for  free  oscillations). 

The  expression  for  a change  in  attenuation  of  a mode  or  a body  wave  is: 


6q  = 


dr  (<5q^M  + Sq^K) 


(IV-17) 


Obviously,  if  the  data  are  available  for  both  Eqs.  (IV-16)  and  (IV-17),  the  two  systems  can 
be  solved  simultaneously.  The  advantage  of  this  procedure  over  solving  for  Eq.  (IV-17)  first 
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and  then  substituting  q (r)  and  qK(r)  into  Eq.  (IV-16)  is  the  limited  resolving  power  of  the 
M’ 


and  the  fact  that  Eq.  (IV-16)  should  impose  additional  constraints  on  q,  if  the  data  cover  rather 
broad  band  of  frequencies  (3-1/2  decades,  in  our  case). 
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Following  the  results  of  Davies  and  Dziewonski,  we  formulate  the  inverse  problem 
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Q data 


for  the 


elastic  parameters  in  terms  of  the  fourth-order  finite-strain  theory.  This  is  a further  develop- 
ment of  a parametric  approach  to  the  description  of  an  earth  model,  first  proposed  by  HAles 

56  5*7 

et  al.  and  implemented  by  Dziewonski  et  al.,  who  derived  a family  of  parametric  eartri  models 

(PEMs).  The  relevant  expressions  are 

p = p0(  l-2€)3/2 


1 2 


= [(1-  Ze)  (L4  + eLz  + \ c 


l3)/p0]1/2 


Vs  = [(1  - 2e)  (M4  4 eM2  4 e2M 3)/p0]l/2 
The  strain  € can  be  represented  as  a function  of  radius: 


<IV-18) 


N 


c(r)  = Yj  anr 
n-0 


2n 


(IV-19) 


By  perturbing  Eqs.  (IV-18)  with  respect  to  e,  zero  pressure  density  pQ,  and  parameters  L 

and  M,  we  express  6p,  SV  , and  6V  in  Eq.  (IV-16)  in  terms  of  the  parameters  of  th< 

P ® 

strain  theory.  A further  assumption  that  a perturbation  in  strain  can  be  represented  by 

N 

2n 


6e(r)  = Yj  bnrZ 
n=0 


(IV-20) 


completes  parametrization  of  the  linearized  equation  (IV-16). 

The  finite-strain  approach  has  been  applied  to  three  regions  of  the  earth's  interior:  lower 
mantle,  outer  core,  and  inner  core.  Because  of  the  superadiabatic  gradients  and  lack  o£  knowl- 
edge of  the  details  of  the  upper-mantle  structure,  the  PEM-like  representation  was  retained 
for  depths  less  than  670  km. 

Table  IV-2  compares  the  strain  function  and  the  finite-strain  parameters  of  the  model 

FSQMK  (finite  strain,  Q , and  QK)  with  those  of  a finite -strain  model  FS,  which  was  constructed 

M" 

ignoring  the  effect  of  velocity  dispersion.  The  density  and  seismic  velocities  in  the  radius  range 
from  0 to  5701  km  can  be  easily  reconstructed  by  substitution  of  appropriate  parameters  into 
Eqs.  (IV-18).  The  seismic  parameters  for  the  upper  mantle  and  the  crust  are  listed  in  Table  IV-3. 
The  attenuation  model  is  presented  in  Table  IV-4;  it  is  quite  similar,  in  principle,  to  that  de- 
scribed by  Sailor  and  Dziewonski  in  Sec.  F above.  For  periods  other  than  1 sec,  the  velocities 
of  model  FSQMK  must  be  corrected  according  to  the  formulas: 


Vg(r,  T)  = Vg(r,  1) 


[ 1 — Q , (r)  • in  T/tt] 


Vp(r,  T) 


Vp(r,  1) 


1 - 


in  T 


4 VS2  -1 

I fz  + 

p 


(1_M)Q;1(r)]l  ' 


(IV-2  l) 


95 


TABLE  IV- 2 

STRAIN  AS  A FUNCTION  OF  SQUARED  RADIUS,  Z = (r/o)2,  AND  FINITE-STRAIN 
PARAMETERS  FOR  TWO  EARTH  MODELS  FSQMK  (FINITE  STRAIN,  Q„,  AND  QK) 
AND  FS  (FINITE  STRAIN,  NO  VELOCITY  DISPERSION)  FOR  THE  INNER  CORE, 
OUTER  CORE,  AND  LOWER  MANTLE.  PARAMETERS  FOR  FSQMK  ARE  VALID 

AT  1-sec  PERIOD. 


Finite-Strain  Parameters 

Region  and 
Radius  Range 
(km) 

Strain 

3 

(g/cm  or  Mb) 

FSQMK 

FS 

FSQMK 

FS 

p0 

6.81054 

6.82124 

h 

2.08429 

2.08069 

-0.27054 

-0.27062 

L2 

-12.11473 

-12.15641 

Inner  Core 
(0  to  1226.6) 

+0.37951  • Z 

+0.38196  • Z 

L3 

6. 99397 

7.07030 

M1 

0. 56027 

0.55099 

m2 

-0.24501 

-0.20576 

M3 

-1.24501 

-1.42100 

p0 

6.51551 

6. 52092 

-0.26742 

-0.26798 

h 

1.26970 

1.22990 

Outer  Core 

(1225.5  to  3484.5) 

+0.31067-  Z 

+0. 30803  • Z 

L2 

-10.  80603 

-10.56040 

+0.15940-  Z2 

+0.16041  • Z2 

L3 

17. 49536 

20.68498 

-0.18517 

-0.18392 

p0 

3. 99624 

3.98614 

+0.21237-  Z 

+0.21180  . Z 

L1 

3.81621 

3.81471 

-0.02661  • Z2 

-0.02802  • Z2 

L2 

-21.85339 

-21.73004 

Lower  Mantle 
(3484.5  to  5701.0) 

L3 

-61.33473 

-59.40508 

M, 

1.24459 

1 . 22405 

m2 

-5.  84662 

-6.02182 

M3 

-36.51793 

-38. 44252 

96 


TABLE  IV-3 

DENSITY  AND  SEISMIC  VELOCITIES  OF  THE  UPPER  MANTLE  AND  THE  CRUST 
OF  EARTH  MODELS  FSQMK  (FINITE  STRAIN,  Qn,  AND  QK;  VALID  AT  1-sec  PERIOD) 

AND  FS  (FINITE  STRAIN,  ATTENUATION  NOT  ACCOUNTED  FOR) 

' L 

Radius 

(km) 

FSQMK 

« i 

Density 

(g/cm3) 

V 

P 

(km/sec) 

V 

s 

(km/sec) 

Density 

(g/cm3) 

V 

P 

(km/sec) 

V 

s 

(km/sec) 

5701 

4.077 

10.034 

5.499 

4.076 

10.009 

5.459 

5801 

3.955 

9. 828 

5.331 

3.954 

9.  806 

5. 290 

5901 

3.830 

9.620 

5.159 

3.  829 

9.601 

5.119 

5951 

3.767 

9.514 

5.073 

3.766 

9. 496 

5.032 

5951 

3.522 

9.198 

4. 908 

3.549 

9.022 

4.818 

6051 

3. 463 

8. 887 

4.  738 

3.489 

8.  703 

4. 645 

6151 

3.403 

8.571 

4.565 

3. 429 

8.378 

4.469 

6151 

3.403 

7.885 

4.400 

3.429 

7.892 

4. 336 

6291 

3.316 

7.885 

4.400 

3. 343 

7.892 

4. 336 

6291 

3.316 

7. 924 

4.  716 

3.343 

7.931 

4. 689 

6352 

3.278 

7.  924 

4.  716 

3.305 

7. 931 

4. 683 

6352 

2.898 

6.500 

3.750 

2.898 

6.500 

3.750 

6357 

2.898 

6.500 

3.750 

2.  898 

6.500 

3.750 

6357 

2.  798 

5.000 

3. 550 

2.798 

6.000 

3.550 

6368 

2.798 

6.000 

3. 550 

2.798 

6.000 

3.550 

6368 

1.030 

1.500 

- 

1.030 

1.500 

- 

6371 

1.030 

1.500 

- 

1.030 

1.500 

- 
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TABLE  IV-4 

Q AND  Qk  IN  THE  EARTH  MODEL  FSQMK. 
Q IS  INFINITE  FOR  DEPTH  INTERVALS  OTHER 
THAN  SPECIFIED  IN  THE  TABLE. 

Depth  Range 

Q 

Q 

(km) 

H 

^ K 

3 to  80 

450 

y 

80  to  420 

93 

412 

420  to  670 

183 

j 

412 

i 

670  to  2886.5 

366 

y 

The  differences  between  FSQMK  and  FS  models  are  similar,  in  general,  to  those  between 

QM1  (Ref.  46)  and  C 2 (Ref.  58).  The  largest  differences  are  in  the  upper  mantle,  particularly  in 

the  depth  range  from  220  to  420  km;  in  our  model,  there  was  also  an  appreciable  increase  in  the 

compressional  velocity  in  this  depth  range.  Overall,  the  increase  in  shear  velocity  in  the  upper 

46 

mantle  is  somewhat  less  than  2 percent  estimated  by  Hart  et  al.  On  the  other  hand,  our  in- 

46 

crease  in  shear  velocity  in  the  lower  mantle  is  roughly  twice  that  of  Hart  et  a^.;  this  is  clearly 
related  to  the  fact  that  our  Q in  this  region  is  only  one-half  of  that  used  to  derive  model  QMl 
(Ref.  46). 

Model  FSQMK  gives  only  marginally  worse  fit  to  the  normal-mode  data  than  other  recent 

earth  models.  The  relative  rms  error  for  1049  normal-mode  data  is  0.20  5 percent  for  FSQMK, 

0.180  percent  for  FS,  0.186  percent  for  PEM-A  (Ref.  57),  and  0.182  percent  for  1066B  (Ref.  28). 

The  error  for  FSQMK  may  be  reduced  in  further  iterations. 

It  is  interesting  to  compare  the  travel  times  computed  for  the  models  FSQMK  and  FS  with 

the  observed  travel-time  data.  In  Fig.  IY-30,  we  show  the  P travel  times  for  a surface  and  a 

deep  focus.  The  baseline  predicted  for  the  model  FSQMK  is  practically  identical  to  that  of  ”1968 

Tables”  (Ref.  47);  in  addition,  for  a focus  at  550  km  depth,  agreement  is  also  reached  with  the 

48 

data  of  Sengupta  and  Julian.  In  Fig.  IY-31,  we  compare  S travel  times;  the  agreement  with 

49 

smoothed  observations  of  Hales  and  Roberts  is  nearly  perfect;  it  is  also  good  for  the  deep- 
focus  data  of  Sengupta,^0  although  there  is  a systematic  difference  of  slightly  over  1 sec  at  dis- 
tances greater  than  7 5°.  As  this  difference  has  the  opposite  sense  from  that  obtained  for  the 
surface  focus,  one  may  suspect  discrepancies  between  the  data  sets.  Also,  the  fit  could  be  im- 
proved by  increasing  the  assumed  period  of  the  deep-focus  S waves. 

What  is  important  is  that  we  have  found  an  excellent  agreement  between  the  baselines  pre- 
dicted from  our  model  and  those  derived  from  the  travel-time  observations  using  land-based 

stations.  We  also  satisfy  within  0.4  sec  the  average  vertical  ScS  travel  time  of  Sipkin  and 
59 

Jordan.  Assuming  that  the  travel-time  baseline  values  are  not  biased,  one  could  conclude  that 
the  differences  between  the  harmonic  average  velocities  under  oceans  and  continents  are  rather 
small;  this  appears  to  agree  with  the  conclusions  of  Okal  and  Anderson. 

A.  M.  Dziewonski 
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I.  CORRELATION  OF  SEISMIC  ACTIVITY  WITH  CHANGES 
IN  THE  RATE  OF  ROTATION  OF  THE  EARTH 

Previous  work  * c suggested  that  seismic  activity  varies  with  time  in  a nonrandorr 

since  it  appears  to  be  possible  to  correlate  these  variations  across  many  of  the  seismic 

of  the  globe.  The  mechanism  for  these  changes  in  seismicity  was  not  discussed  in  these 

though  it  is  clear  that  only  global  mechanisms  are  candidates,  and  from  the  rather  limited  num- 

6 3 

ber  of  possibilities  the  rate  of  rotation  of  the  earth  must  be  a leading  contender. 

Detailed  analysis  of  earthquake  time-series  from  many  different  areas  indicated  th4t  the 
following  features  appear  quite  commonly:  A small  trough  in  seismicity  during  late  1964  is  fol- 
lowed by  a high  during  1965,  and  then  a period  of  unusually  low  activity  extends  from  mi|d-19  66 
through  most  of  1967.  Activity  commonly  increases  during  19  68,  and  then  becomes  mudh  less 
coherent  from  the  period  1969  to  1973,  though  there  is  a suggestion  of  an  increase  in  activity 
in  1971. 

Changes  in  the  rate  of  rotation  of  the  earth  (i.e.,  in  the  length  of  the  day)  have  been  (well- 
documented  since  19  55,  when  atomic  clocks  were  introduced.  Figure  IY-32  shows  the  length  of 


the  day  from  mid-1955  to  mid-1976,  using  data  sampled  at  30-day  intervals,  supplied  to 


us  by 


Dr.  L.  V.  Morrison  (Royal  Greenwich  Observatory,  private  communication).  A very  strfong  sea- 
sonal term  with  amplitude  approximately  1 msec  and  containing  6-  and  12 -month  components^* 
is  apparent.  These  seasonal  terms  are  effectively  removed  by  passing  a sliding  12-month  aver- 
aging window  through  the  data.  The  resulting  curve  is  labeled  " Annual  Means"  in  Fig.  IY-32. 

This  mean  variation  in  rotation  rate  has  a great  amount  of  structure.  In  particular*  notice 
the  rapid  increases  in  the  length  of  day  (corresponding  to  decelerations  of  the  earth)  during 

1956-57,  1963-64,  and  1971.  Smaller  increases  are  clear  in  1965  and  1968. 

6 3 

Anderson  0 has  suggested  that  the  large  number  of  earthquakes  in  the  period  1900  tb  1910 
may  have  been  associated  with  the  very  rapid  change  in  the  length  of  day  that  began  aboUt  1890. 

It  therefore  seems  logical  to  compare  seismicity  with  the  rate  of  change  of  the  length  of  day  in 
Fig.  IV-32.  Interestingly,  the  periods  of  most-rapid  change  all  seem  to  be  associated  with  high 
seismic  activity.  Figure  IV-3  3 shows  counts  of  large  earthquakes  during  the  period  19  95  to  1976, 
and  peaks  in  19  56-57  and  1971  are  noticeable.  There  is  no  clear  peak  associated  with  the  1963-64 
change,  though  this  culminated  in  two  large  earthquakes  — the  Kurile  Island  (1963)  and  Alaska 
(1964)  events.  Figure  IY-34  shows  seismic-energy  release  during  the  same  period,  and  now  all 
three  peaks  emerge.  However,  the  calculation  of  energy  requires  reliable  estimates  of  the  mag- 
nitude of  large  events,  and  this  is  notoriously  difficult.  In  particular,  it  has  been  necessary  to 
6 5 

use  the  Duda  catalog  for  the  period  19  55  to  1964,  and  it  is  not  clear  that  the  listed  magnitudes 
in  this  catalog  are  very  reliable. 

Comparing  the  rate  of  change  of  the  length  of  day  with  the  earthquake  time  series  described 
earlier  shows  an  interesting  agreement.  In  particular,  the  change  in  rotation  rate  was  {small 
during  the  latter  part  of  1964,  and  was  very  low  during  the  period  1966-67,  as  was  the  Seismicity. 

The  preliminary  conclusion  of  this  study  is  that  there  seems  to  be  sufficient  correlation  be- 
tween changes  in  seismic  activity  and  changes  in  the  rate  of  rotation  of  the  earth  to  warrant  fur- 
ther investigation.  It  is  interesting  to  note  that  the  rate  of  rotation  is 'decreasing  during  1976, 
perhaps  at  a significantly  fast  rate.  We  also  note  that  there  have  been  a number  of  large  earth- 
quakes during  1976.  This  somewhat  speculative  correlation  will  be  examined  further  as  more 

rotation-rate  data  become  available.  A 

M.  A.  Chinnery 
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COHERENCE 


Fig.  IV-1.  Coherence  spectra  of  correlated 
white  processes.  B3  is  Burg  Estimator  of 
order  3;  T3  is  Taper  and  Transform  Esti- 
mator of  lag  3;  T50  is  Taper  and  Transform 
Estimator  of  lag  50. 


Fig.  IV-3.  Cross-spectral  amplitudes 
of  single  sinusoid.  B4  is  Burg  Estima- 
tor of  order  4;  T4  is  Taper  and  Trans- 
form Estimator  of  lag  4;  T50  is  Taper 
and  Transform  Estimator  of  lag  50. 


Fig.  rV-2.  Phase  spectra  of  correlated 
white  processes.  B3  is  Burg  Estimator 
of  order  3;  T3  is  Taper  and  Transform 
Estimator  of  lag  3;  T50  is  Taper  ,and 
Transform  Estimator  of  lag  50. 


Transform  Estimator  of  lag  50. 
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COHERENCE 


Fig.  IV- 5.  Coherence  spectra  generated  by  autoregressive  estimators  of  orders  2, 
5,  10,  15,  20,  30,  40,  50. 
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Fig.  IV-6.  Phase  spectra  generated  by  autoregressive  estimators  of  orders  2 
5,  10,  15,  20,  30,  40,  50. 


105 


FREQUENCY  (proportion  of  Nyquist) 


FREQUENCY  (proportion  of  Nyquist) 


Fig.  IV-7.  Coherence  and  phase  spectra  for  varying  SNR. 
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Fig.  IV-8.  20  March  1976  event  recorded  at  NORSAR,  GBA,  YKC,  and  COL. 
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Fig.  IV-9.  Sum  of  eight  identical  SP  waveforms  to  simulate  pP. 
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Fig.  IV-11.  Surface  waves  recorded  at  SRO,  MAIO  for  20  March  1976  event 
and  4 July  1976  explosion  with  same  epicenter. 


Fig.  IV -12.  Mg-mt,  for  20  March  1976  event 
(point  A)  superimposed  on  Marshall  and 
Basham's17  Mg-m^  data  for  Kazakhstan. 
Assuming  4 multiple  events  contribute  to 
surfaces  equally,  Ms-mb  for  each  event 
moves  to  point  B.  Assuming  average  body- 
wave  amplitude  per  event  is  1.5  bigger  than 
first  event,  Ms-m^  falls  squarely  on  ex- 
plosion trend  line  (point  C). 
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Fig.  IV-13.  First  4 sec  of  SP  seismogram 
recorded  at  Yellowknife  for  20  March  1976 
event.  Arrows  indicate  interpretation  of 
cepstrum. 
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Fig.  IV-14.  Cepstrum  of  seismogram 
in  Fig.  IV-13  indicating  two  events,  the 
second  larger  than  the  first  and  arriv- 
ing about  0.9  sec  later.  Both  depth 
phases  are  at  approximately  0.8  sec, 
indicating  very  shallow  (a  kilometer  or 
so)  events. 
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Fig.  IV- 15.  Kurile  Island  event. 
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Fig.  IV-16.  Low-pass  cepstrum 
of  Kurile  Island  event. 


Fig.  IV- 17.  Cepstrum  of  Kurile  Islan4  event. 
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Fig.  IV-18.  MEM  cepstrum  of  order  60. 


Fig.  IV- 19. 
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Fig.  IV-20.  ME  cepstrum 
of  pP  portion  only. 
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Fig.  IV-21.  Window  of  3-component  data  and  its  corresponding  particle -motion 
ellipsoid  which  is  least-squares  approximation  to  actual  particle  motion.  Ellipsoid 
axes  have  lengths  if,  i 2,  and  i3,  where  I * is  called  principal  component  of  data. 
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Fig.  rV-22.  Example  of  polarization  filtering  applied  to  an  earthquake  in  Kuriles,  mb  5.4,  recorded  at  MAIO  66°  away. 
Traces  1 to  3 are  Z,  R,  and  T components  of  unfiltered  data  beginning  at  origin  time.  Traces  4 to  6 are  same  compo- 
nents filtered  to  pass  periods  from  11.11  to  33.33  sec.  Trace  7 is  instantaneous  gain  factor  used  in  polarization  scheme. 
Traces  8 to  10  are  polarized  traces  obtained  by  projecting  principal  component  of  ellipsoid  onto  data,  scaled  by  gain  factor. 


Fig.  IV-23.  Two-component  particle  motions  for  (a)  idealized  P waves 
and  (b)  Rayleigh  wave.  If  R component  of  surface  wave  is  Hilbert  trans- 
formed, particle  motion  is  changed  from  elliptically  to  linearly  polarized. 
For  P waves,  motion  is  converted  from  linear  to  elliptical. 
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Fig,  IV-24.  Example  of  polarization  filtering  applied  to  two  components  Z and  R of  Fig.  IV-22.  Traces  1 and  2 are 
bandpassed  components.  Trace  4 is  R component  after  being  Hilbert  transformed.  Using  a gain  based  on  ellipticity 
only,  trace  5,  surface  waves  can  be  passed  or  rejected,  as  shown  respectively  by  trace  pairs. 6 and  7,  8 and  9.  An 
improyed  gain  factor  given  by  correlation  coefficient,  trace  10,  produced  better  rejection  of  body  waves  in  traces  11 
and  12,  while  passing  surface-wave  motion.  Traces  13  and  14  preserve  body  phases  while  surface  waves  are  sub- 
tracted from  data. 


Fig.  IV-25.  Attenuation  of  0S0  following  1964  Alaskan  earthquake  as  observed 
on  UCLA  gravimeters  #4  and  #7.  Each  point  represents  power  in  decibels  in 
a different  6 -day  window,  centered  at  the  value  given  on  horizontal  axis.  Max- 
imum value  of  power  in  first  window  (0  to  6 days)  is  given  reference  value  of 
100  dB.  Range  of  Q given  corresponds  to  plus-or-minus  one  standard  devia- 
tion of  l/Q  as  obtained  from  regression. 


Fig.  IV-26.  Attenuation  of  0S2  singlets 
and  integrated  multiplet  on  UCLA  gra- 
vimeter #4.  Q can  be  measured  sep- 

_ A 

arately  for  two  distinct  peaks  0S2  and 
+ 1 

0S2  . For  explanation,  see  text. 
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Fig.  IV-27.  Two  models  of  mantle,  MM8  and  RW1.  Model  MM8  of 
Anderson  et  al.33  Was  used  as  a starting  model  in  "ranking  and  winnow- 
ing" data-space  inversion  to  obtain  model  RW1.  We  used  observed 
values  and  corresponding  standard  errors  for  a representative  set  of 
40  modes  (3  toroidal,  7 radial,  aSj,  8 high-Q  overtones,  and  21  qS^). 
Relative  rms  residual  for  model  MM8  is  0.416;  this  was  reduced  to 
0.196  for  final  model,  RW1. 
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Fig.  IV-28.  Comparison  of  velocity 
dispersion  associated  with  a "trian- 
gular” attenuation  function  and  that 
due  to  an  absorption-dispersion  pair.41 


Fig.IV-29.  Velocity  dispersion  due 
to  four  attenuation  functions.  Note 
that  dispersion  curves  correspond- 
ing to  attenuation  functions  la  and  lb 
as  well  as  2a  and  2b  are  very  simi- 
lar, respectively. 
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Fig.  IV- 30.  Comparison  of  observed  and 
computed  travel  times  for  P waves. 
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Fig.  IV- 31.  Comparison  of  observed  and 
computed  travel  times  for  S waves. 
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Fig.  IV-32.  Length  of  day  from  mid-1965  to  mid-1976,  plotted 
at  30-day  intervals. 
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Fig.  IV-33.  Comparison  of  rate  of  change  of  length  of  day  with  counts 
of  large  earthquakes. 
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Fig.  IV-34.  Comparison  of  rate  of  change  of  length  of  day  with  seismic -energy 
release  in  6 -month  intervals. 
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V.  DATA  SYSTEMS 


A.  PDP-11  COMPUTER  FACILITY 

The  PDP-11  computer  is  presently  configured  as  a time-sharing  system  running  the  UNIX 
operating  system.  The  central  processor  is  a PDP-11 -50  with  112K  words  of  memory.  Also 
included  in  the  configuration  are  a number  of  terminals,  on-line  and  off-line  data-storagd  facil- 
ities, a Versatec  printer  plotter,  a card  reader,  a floating-point  processor,  and  a connection  to 
the  ARPANET.  Terminals  include  4 Tektronix  4014  graphics  terminals,  5 DEC  VT52  Yideo 
terminals,  5 DEC  LA36  hard-copy  terminals,  and  1 Diablo  typewriter  terminal.  Data-storage 
’devices  include  4 disk  drives  and  2 tape  drives.  There  are  two  types  of  disk  drives:  two 
40M-byte  DIVA  DD52  drives  used  for  user  file  storage,  and  two  2.5M-byte  DEC  RK05  drives 
used  for  system  storage  and  swapping.  The  tape  drives  include  a 7-  and  a 9 -track  drive.  This 
system  supports  a user  community  of  about  40  people  doing  a variety  of  computer  tasks.  These 
tasks  include  some  limited  Fortran  number -crunching  programs,  computer  graphics  in  bbth 
Fortran  and  C,  access  to  the  ARPANET,  and  the  development  of  user  and  system  software. 

In  addition  to  normal  software  maintenance  and  the  incorporation  of  externally  acquired 
software  into  our  system,  there  has  been  significant  progress  on  our  graphics  software,  )pur 
software  to  access  the  datacomputer,  and  our  methods  of  locally  dealing  with  SRO  data.  ^The 
following  information  is  concerned  with  all  these  areas. 

The  UNIX  graphics  system  is  a device-independent  user-oriented  graphics  system.  It  is 
designed  to  be  implemented  for  a variety  of  devices  ranging  from  Tektronix  displays,  to  Refresh 
scopes,  to  Versatec  plotters.  The  device  independence  is  accomplished  through  a universal 
graphics  language  (UGL)  which  is  based  on  the  network  graphics  protocol  as  described  in  the 
ARPANET  Protocol  Handbook.  UGL  allows  the  user  to  execute  a set  of  graphics  primitives 
which  combine  to  form  a display  structure.  The  structure  represents  a graphics  picture  which 
can  be  manipulated  with  additional  graphics  primitives.  The  graphics  primitives  are  imple- 
mented through  a set  of  subroutines  for  both  Fortran  and  C. 

The  basic  conceptual  building  block  of  the  UNIX  graphics  system  is  the  picture  segment.  A 
picture  segment  is  a logically  related  grouping  of  vectors,  points,  and  character  strings  Kvhich 
represent  a portion  of  the  picture.  A complete  picture  can  therefore  be  thought  of  as  one  or 
more  picture  segments. 

Manipulation  of  a picture  is  accomplished  through  the  execution  of  a number  of  graphics 
primitives  on  individual  picture  segments.  Segments  may  be  created,  destroyed,  duplicated, 
translated,  or  scaled  independently  in  x and  y.  They  may  also  be  displayed  (posted)  or  kept 
hidden  (unposted)  in  the  actual  picture  being  displayed.  An  unposted  segment  is  still  logically 
part  of  the  picture,  although  it  is  not  actually  visible. 

The  manipulation  of  picture  segments  and  the  retention  of  the  entire  picture  in  a data  file 
allow  the  user  to  change  portions  of  a picture  without  recomputing  the  unchanged  portions.  In 
addition,  the  stored -picture  data  file  also  allows  the  user  to  move  a created  picture  from  one 
graphics  device  to  another  by  simply  moving  the  data  file.  A special  command  is  provided 
which  reads  back  the  existing  graphics  data  file  from  one  device  and  sends  it  to  another  graphics 
device. 

All  the  overhead  of  maintaining  the  picture  data  file  and  the  manipulation  of  segment^  is 
done  by  a process  separate  from  the  user  program.  This  frees  the  user  from  having  to  allocate 
and  maintain  storage  space  to  support  the  graphics-picture  file. 
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In  some  cases,  there  are  special  features  associated  with  a graphics  device  which  are  not 
implemented  in  the  hardware  of  the  other  devices,  e.g.,  the  multiple  character  sizes  and  the  ma- 
nipulation of  a cursor  on  the  Tektronix.  We  decided  that  these  features  should  not  be  sacrificed 
for  the  sake  of  compatibility.  The  commands  which  implement  these  features  are  either  simu- 
lated in  software  for  the  other  devices  or  result  in  a warning  message  to  the  user's  error  file. 

Graphics  is  implemented  through  the  use  of  graphics  daemons  (free -running  UNIX  processes) 
which  are  always  present  in  the  system.  Each  graphics  device  has  its  own  daemon  which  is  re- 
sponsible for  maintaining  the  software  overhead  for  that  device.  The  user  communicates  with  a 
graphics  daemon  through  a set  of  subroutines  which  send  graphics  commands  to  the  desired 
graphics  daemon  using  UNIX  interprocess  as  communication-channels  called  ports.  Each  daemon 
sets  up  an  input  and  an  output  port  during  system  initialization  and  makes  them  accessible  to  all 
users  in  the  system.  Each  graphics  subroutine  sends  a message  over  the  port  corresponding  to 
a graphics  command  and  its  parameters.  The  graphics  daemon  receives  the  command  messages 
and  stores  them  in  a display  file.  The  display  file  is  a structured  file  which,  in  addition  to  the 
graphics  commands,  maintains  the  information  necessary  to  identify  picture  segments  and  their 
status.  This  is  done  with  a picture-segment  table  containing  pointers  to  buffers  of  device- 
dependent graphics  commands  and  with  flags  indicating  the  status  of  each  segment.  The  display 
file  also  contains  the  x and  y translation  and  scaling  for  each  picture  segment.  The  paint 
command  causes  the  graphics  daemon  to  transform  the  display  file  into  device -dependent  com- 
mands, and  to  transmit  these  commands  to  the  appropriate  graphics  device. 

The  currently  implemented  graphics  system  is  limited  to  Tektronix  graphics.  Work  on  a 
Versatec  plotter  is  about  70 -percent  complete.  There  are  no  plans  to  implement  the  graphics 
system  for  any  other  specific  devices  at  the  present  time,  although  the  design  does  allow  for 
future  expansion. 

Our  basic  access  to  the  datacomputer  is  presently  by  means  of  DCI,  which  is  a new  program 
that  provides  a convenient  terminal  interface  between  our  local  UNIX  system  and  the  datacom- 
puter. This  program  was  designed  to  work  in  either  of  two  basic  situations: 

(1)  The  user  has  a prepared  datalanguage  "program"  which  he  wishes  to 
transmit  to  the  datacomputer  with  a minimum  of  effort. 

(2)  The  user  wishes  to  work  with  the  datacomputer  in  an  interactive  manner, 
such  as  when  he  is  trying  out  new  sequences  of  datalanguage. 

Some  of  the  main  features  of  DCI  are  summarized  in  the  following  paragraphs. 

The  datacomputer  sends  the  user  a great  variety  of  messages.  Because  of  the  large  volume 
of  these  messages,  DCI  allows  the  user  to  choose  between  three  different  verbosity  levels:  un- 
der FULL,  all  messages  will  be  typed  out;  under  NORMAL,  only  error  messages  and  urgent 
informational  messages  will  be  output;  and  under  SILENT,  all  messages  will  be  suppressed. 

When  the  user  has  a sequence  of  datalanguage  and  DCI  commands  that  will  be  used  frequently, 
he  may  create  a UNIX  file  that  contains  these  commands,  and  then  direct  DCI  to  take  its  input 
from  this  file.  A command  file  may  also  contain  comment  lines  and  references  to  other  com- 
mand files  up  to  three  levels  of  nesting.  Another  verbosity  level,  called  MONITOR,  will  cause 
each  line  from  the  command  file  to  be  typed  out  on  the  user's  terminal  as  it  is  processed  by  DCI. 

Using  DCI,  data  maybe  transferred  in  either  direction  between  the  datacomputer  and  a UNIX 
file,  using  either  the  default  datalanguage  connection  or  a separate  data  connection.  The  user 
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may  establish  a separate  connection  by  specifying  only  the  direction  of  transfer,  the  byte  size, 
and  the  datacomputer  port  involved.  DCI  handles  the  details  of  selecting  a socket  number  and 
informing  the  datacomputer  of  the  connection. 

DCI  will  write  all  messages  and  commands,  regardless  of  verbosity  level,  to  a file  galled 
dci. script.  This  will  provide  the  user  with  a complete  record  of  his  datacomputer  sessiqn. 

This  record  is  easily  discarded  if  it  is  not  needed  by  the  user,  but  is  an  invaluable  sourcb  of 
information  when  there  are  problems  in  a transaction. 

While  DCI  is  running,  the  user  may  request  status  information  which  will  include  whjat  state 
the  program  is  in,  what  separate  connections  exist  (if  any),  and,  if  data  are  being  transferred, 
the  elapsed  time  of  the  transfer  and  the  number  of  bytes  that  have  been  passed.  The  use*  may 
also  return  to  the  UNIX  shell  and  edit  command  files  or  inspect  the  script  file  before  returning 
to  DCI. 

The  third  item  to  discuss  is  our  local  handling  of  SRO  data.  Much  of  our  SRO  data  a^d  some 
other  seismic  data  are  currently  on  magnetic  tape  in  waveform  data  unit  format.  The  format 
assigns  a tape  file  to  each  seismogram.  The  file  contains  the  seismogram  plus  alphanumeric 
information  describing  the  seismogram.  As  discussed  in  our  previous  SATS,  this  tape  format 
is  quite  general  and  is  the  primary  medium  for  the  transfer  of  waveform  data  between  oufr  PDP-7 
and  PDP-11  systems.  Two  PDP-11  programs  have  now  been  written  to  exchange  data  on  such 
tapes  with  a waveform  database  which  can  be  kept  in  the  UNIX  system.  Within  UNIX,  thq  in- 
dividual seismograms  continue  to  reside  in  separate  files,  but  the  alphanumeric  information  is 
combined  into  a single  index  file  on  the  disk.  That  file  contains  the  identification  information 
from  the  waveform  data  unit  tape  plus  the  UNIX  file  names  of  the  files  containing  the  adtual 
seismograms.  The  seismograms  themselves  are  stored  as  binary  floating-point  numberls. 

The  two  programs  which  have  been  written  are  DUREAD  and  DUWRITE.  DUREAD  reads 
all  or  part  of  a waveform  data  unit  tape  and  either  adds  the  data  to  an  existing  database  or  will 
create  a new  one.  On  input,  data  are  converted  from  binary-gain-ranged  or  18-bit  binary  inte- 
gers into  PDP-11  floating  point.  The  DUWRITE  program  outputs  all  of  a specified  database 
onto  a waveform  data  unit  tape.  All  alphanumeric  information  will  be  preserved  but  the  repre- 
sentation of  data  samples  is  modified  from  floating-point  to  16-bit  integer  format  on  the  tape, 
so  the  user  must  take  care  to  avoid  having  sample  values  which  are  too  large.  These  programs 
are  temporary  until  more  sophisticated  ones  can  be  developed. 

L.  J.  Turek 
J.  Sax 


INTERACTIONS  WITH  THE  DATACOMPUTER 

We  have  started  to  use  our  datacomputer  interface  program  DCI  to  store  data  on  the 


data- 


computer and  to  fetch  data  from  it.  Three  current  applications  of  DCI  are  described  in  tins 
section.  Other  applications  are  being  developed,  as  are  alternative  software  modules  to  fetch 
and  store  data. 

The  ISC  seismic  bibliography  for  the  years  1965-1973  has  been  stored  on  the  datacoiinputer 
under  the  seismic. lincoln.biblio  node.  This  database  contains  information  about  14,230  pub- 
lished seismic  works:  accession  number,  title,  joumal/publisher  code,  volume,  numbeh,  year, 
page,  authors,  and  keywords.  The  file  is  inverted  on  keyword  value  and  author  surname  so  that 
searches  involving  these  items  will  be  optimized.  A local  UNIX  program,  LOOKUP,  will  use 
DCI  to  connect  to  the  datacomputer,  log  into  the  biblio  node,  and  open  the  appropriate  fil^s  and 
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output  ports.  When  initialization  is  complete,  it  will  accept  any  user  lookup  requests  that  can 
be  expressed  in  datalanguage.  When  the  results  are  received,  the  output  may  be  typed  out  on 
the  user's  terminal  or  written  into  a UNIX  file. 

There  now  exists  in  the  datacomputer  a preliminary  event  summary  file  (PESF)  containing 
event  locations  obtained  by  using  the  large  arrays,  and  selected  SP  stations.  The  PESF  file  is 
updated  frequently  and  is  about  7 to  8 weeks  behind  real  time  for  bulletin  listings.  These  pre- 
liminary epicenters  are  very  useful  to  us  for  determining  the  time  periods  of  interest  for  the 
SRO  digital  data.  The  PESF  file  is  created  and  maintained  by  the  Seismic  Data  Analysis  Center 
and  is  suitably  protected  against  unauthorized  access.  On  a routine  basis,  we  read  data  from 
the  PESF  files  and  transfer  the  data  into  our  UNIX  system.  There  is  a small  amount  of  re- 
formatting that  occurs  during  transfer  between  the  mass  store  and  UNIX.  Once  the  data  are 
within  UNIX,  they  are  reformatted  again  and  additional  information  such  as  geographic  region 
name  and  seismic  region  name  is  added  to  the  list.  The  files  are  then  divided  into  monthly 
periods  and  maintained  on-line  for  immediate  local  access. 

The  International  Seismological  Center  (ISC)  also  makes  available  computer  tapes  containing 
event  lists  and  the  individual  reports  from  stations  for  these  events.  This  is  a large  database 
which  is  difficult  to  use  because  of  its  size  and  the  complexities  and  variabilities  of  its  formats. 
We  are  currently  in  the  process  of  storing  these  data  on  the  datacomputer  for  the  time  period 
1964  to  1973.  The  event  information  is  already  entered,  and  the  supporting  station  data  are  now 
being  entered.  When  these  data  are  entered,  they  can  be  of  use  to  all  members  of  our  group 
with  greater  convenience  than  would  be  possible  were  the  data  kept  on  tape  or  in  some  other  com- 
puter system  with  a large  amount  of  disk  storage. 

The  datalanguage  file  description  for  the  ISC  epicenter  database  has  several  inverted  con- 
tainers. Some  of  the  inverted  containers  are  year,  magnitude,  hemisphere,  and  geographic 
region  numbers  along  with  a few  other  pointers,  including  a sequence -number  counter.  The 
corresponding  event  station  readings  are  loaded  into  a second  database,  one  file  per  year.  Since 
each  event  has  a sequence  number,  this  number  is  carried  over  into  the  station-reading  data- 
base. This  event  number  is  the  pointer  from  one  file  into  the  other.  The  ten  station-readings 
files  have  identical  format  descriptions,  thus  enabling  a user  to  treat  them  as  a group  (a  special 
datacomputer  feature)  when  performing  searches  on  this  database.  The  volume  of  data  con- 
tained within  the  station-readings  file  is  far  greater  than  the  data  contained  in  the  epicenter 
file.  The  average  number  of  readings  per  year  in  the  station -readings  file  is  over  400,000. 

The  total  volume  of  bits  stored  for  the  station  readings  will  eventually  exceed  1400  megabits 
for  the  10 -year  database.  Some  of  the  containers  that  are  inverted  in  the  station-readings  list 
are  station  name,  phase  code,  and  event  number.  This  event  number  is  the  pointer  back  into 
the  event  database  for  the  correct  event  parameters  that  are  associated  with  the  station  reading. 

R.  M.  Sheppard 

L. J.  Turek 
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